AD 


Award  Number:  W81XWH-11-1-0493 


TITLE:  Development  of  Magnetic  Resonance  Imaging  Biomarkers  for 
Traumatic  Brain  Injury 

PRINCIPAL  INVESTIGATOR:  E  Mark  Haacke 


CONTRACTING  ORGANIZATION: 

Wayne  State  University,  Detroit,  Michigan  48201 

REPORT  DATE: 

July  2013 

TYPE  OF  REPORT:  Final 


PREPARED  FOR:  U.S.  Army  Medical  Research  and  Materiel  Command 
Fort  Detrick,  Maryland  21702-5012 


DISTRIBUTION  STATEMENT:  (Check  one) 

XD  Approved  for  public  release;  distribution  unlimited 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 
data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202- 
4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1 .  REPORT  DATE  2.  REPORT  TYPE 

July  201 3  Final 

3.  DATES  COVERED 

2  June  2011-1  June  2013 

4.  TITLE  AND  SUBTITLE 

Development  of  Magnetic  Resonance  Imaging  Biomarkers  for  Traumatic  Brain  Injury 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

W81XWH-1 1-1 -0493 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

E  Mark  Haacke 

email:  nmrimaging@aol.com 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Wayne  State  University,  Detroit,  Michigan  48201 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  Medical  Research  and  Materiel  Command 

Fort  Detrick,  Maryland  21702-5012 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  Public  Release;  Distribution  Unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

The  purpose  of  this  research  program  is  to  develop  an  imaging-based  protocol  to  improve  diagnosis  and  outcome  prediction  of 
mild  traumatic  brain  injury  (mTBI).  Our  major  findings  over  the  last  year  include:  First,  we  have  demonstrated  that  these  new 
advanced  MRI  methods  nicely  complement  conventional  imaging  methods  in  their  ability  to  detect  mTBI  in  the  acute  setting. 
Second,  susceptibility  weighted  imaging  and  mapping  (SWIM)  are  able  to  delineate  venous  structures  and  microbleeds  making 
it  possible  to  quantify  and  monitor  evolution  of  pathological  changes.  Third,  we  have  continued  to  use  DTI  data  to  look  for  local 
variations  in  ADC  and  FA  values  to  assess  tissue  damage.  And  fourth,  we  have  begun  to  accelerate  the  data  collect  process  in 
patients  with  added  support  in  both  Radiology  and  Emergency  medicine  for  both  volunteers  and  mTBI  patients.  The 
significance  of  this  work  to  date  is  that  we  are  now  able  to  find  evidence  of  damage  in  mTBI  where  that  has  heretofore  been 
difficult  with  conventional  MRI.  Keywords:  Mild  traumatic  brain  injury,  susceptibility  weighted  imaging  and  mapping,  temporal 
nature  of  TBI  recovery,  medullary  vein  damage 

15.  subject  terms-  none  provided 

16.  SECURITY  CLASSIFICATION  OF: 


b.  ABSTRACT 


17.  LIMITATION 
OF  ABSTRACT 


18.  NUMBER 
OF  PAGES 


19a.  NAME  OF  RESPONSIBLE  PERSON 

USAMRMC 

19b.  TELEPHONE  NUMBER  (include  area 


U 


code) 


a.  REPORT 

U 


c.  THIS  PAGE 

u 


uu 


96 


Table  of  Contents 


Page 


Introduction .  4 

Body .  5 

Key  Research  Accomplishments .  13 

Reportable  Outcomes . 13 

Conclusion . 14 

Appendices . Attachment 


The  views,  opinions  and/or  findings  contained  in  this  report  are 
those  of  the  author (s)  and  should  not  be  construed  as  an  official 
Department  of  the  Army  position,  policy  or  decision  unless  so 
designated  by  other  documentation. 


INTRODUCTION 


Background:  Traumatic  brain  injury  (TBI)  has  been  labeled  as  a  “signature  wound”  in  the  anti-terrorism 
wars  in  Iraq  and  Afghanistan.  15%  -  25%  of  surveyed  returning  service  members  have  been  reported  to 
have  possible  long-term  mild  traumatic  brain  injury  (mTBI)  or  concussion.  In  the  civilian  sector,  the 
prolonged  neuro-cognitive  and  functional  symptoms  following  mTBI  affects  over  1.2  million  Americans 
annually.  However,  our  understanding  of  the  neuropathology  of  mTBI  and  its  recovery  process  is  still 
very  limited  due  to  the  lack  of  sensitive  clinical  and  diagnostic  tools.  First,  current  TBI  classification 
scales  are  mainly  based  on  results  from  computed  tomography  (CT)  scans  of  moderate  to  severe  TBI 
patients,  and  thus  can  hardly  be  applied  to  mTBI  cases.  Secondly,  current  laboratory  biological  markers, 
clinical  CT  and  conventional  magnetic  resonance  imaging  (MRI)  are  either  insensitive  or  not  specific  to 
the  subtle  abnormalities  in  mTBI  and  poorly  predict  patient’s  long-term  outcome.  It  is  in  urgent  need  of 
developing  a  set  of  advanced  MR  imaging  biomarkers  that  can:  i)  be  sensitive  enough  to  differentiate 
mTBI  patient  population  from  normal  healthy  subjects,  ii)  have  the  potential  for  outcome  prediction;  and 
iii)  assist  in  management  of  mTBI  patients  in  acute  settings. 

Objective/Hypothesis:  The  overall  goal  of  the  current  research  is  to  develop  imaging-based  biomarkers 
to  improve  diagnosis  and  outcome  prediction  of  mTBI.  Advanced  MRI  techniques  reveal  many  more 
details  than  conventional  imaging.  For  example,  susceptibility  weighted  imaging  (SWI)  and  susceptibility 
mapping  (SWIM)  can  detect  and  quantify  temporal  reduction  of  hemorrhagic  lesions  associated  with 
mTBI  patients’  functional  outcome  and  diffusion  tensor  imaging  (DTI)  can  detect  axonal  injury  at  the 
acute  stage  and  possible  changes  over  time.  As  the  next  generation  of  SWI  developed  in  our  lab,  SWIM 
will  be  used  to  quantify  iron  in  microbleeds  and  oxygen  saturation  in  major  veins  throughout  the  brain. 
Our  central  hypothesis  is  that  axonal  injury  (measured  by  DTI)  and  vascular  damage  (detected  as 
hemorrhagic  lesions  by  SWI/SWIM)  are  important  pathologies  in  mTBI  that  are  associated  with  patients’ 
neurocognitive  and  clinical  symptoms  in  their  recovery. 

Study  Design:  Adult  mTBI  patients  are  screened  and  enrolled  at  the  acute  stage  (within  24  hours  after 
injury)  from  the  emergency  department  of  our  local  Level-One  Trauma  Center.  Both  up-to-60-minute 
MRI  scan  and  neuropsychological/clinical  assessment  will  be  conducted  in  mTBI  patients.  In  the  subacute 
(1  month)  and  chronic  (6  months)  stages  after  injury,  the  patients  arebrought  back  to  repeat  both  the  MRI 
scans  and  neurocognitive  evaluations.  Age/gender/education-matched  healthy  controls  will  be  recruited 
and  followed  up  using  the  same  imaging  and  neuropsychological  protocol  of  evaluation. 

Specific  Aims:  Specific  Aim  1  is  to  assess  whether  the  advanced  MRI  data  (SWI  and  DTI)  acquired 
acutely  (within  24  hours  after  injury)  are  more  sensitive  than  conventional  imaging  (CT  and  cMRI)  in 
detecting  mTBI.  Specific  Aim  2  is  to  determine  at  what  time  point  after  the  injury  (i.e.,  acute,  subacute 
and  chronic)  these  advanced  MR  techniques  can  differentiate  mTBI  patients  from  healthy  controls. 
Specific  Aim  3  is  to  use  susceptibility  weighted  imaging  and  mapping  (SWIM),  our  next  generation  of 
SWI  technique,  to  quantify  the  amount  of  iron  in  microbleeds,  to  monitor  any  changes  (evolution)  of  the 
microbleed  over  time  and  to  monitor  oxygen  saturation  compared  to  normal  volunteers.  And  Specific 
Aim  4  is  to  determine  which  MRI  indices  have  statistically  significant  associations  with  neurocognitive 
outcome  in  mTBI  patients  over  all  time  points. 
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BODY 


Summary  of  the  Statement  of  Work  for  Phase  1 : 

Phase  I  (Years  1-2):  Recruit  a  cohort  of  50  mTBI  patients  from  the  emergency  setting  of  our  local  Level- 
One  trauma  center,  and  50  non-TBI  volunteers  to  establish  and  validate  an  image-based  diagnostic 
approach  to  brain  injury. 

a)  Screen  and  enroll  mTBI  patients;  perform  the  MRI  protocol  (conventional  imaging,  DTI  and  SWI); 
conduct  NP  tests  and  clinical  outcome  measures. 

Due  to  the  late  approval  of  IRB,  we  started  enrollment  fairly  late  (March,  2012).  To  date,  we  have 
successfully  recruited  26  mTBI  patients  and  scanned  them  at  the  acute  stage  (within  24  hours  after 
injury).  Our  momentum  dramatically  improved  this  last  quarter  and  we  have  strong  momentum  now  that 
we  are  confident  will  make  it  possible  to  easily  complete  the  study  over  the  next  year. 

b)  Screen  and  enroll  healthy  volunteers  and  perform  the  same  procedures  as  mentioned  above. 

We  have  successfully  enrolled  22  control  subjects  in  the  past  year.  As  above,  our  momentum  here  has 
dramatically  improved  and  we  are  getting  key  data  here  for  normal  controls. 

c)  Follow  up  mTBI  patients  at  1 -month  and  6-months  after  injury  and  controls  at  1 -month  after  first  scan 
by  using  the  same  imaging  protocol,  NP  tests  and  outcome  measures. 

We  have  successfully  performed  1 -month  follow  up  in  our  21  patients  and  21  controls.  We  also  followed 
up  in  all  4  patients  who  are  mature  for  6-month  follow  up. 

d)  Analyze  data  and  establish  a  quantitative  approach  to  identifying  abnormal  brain  voxels  on  DTI  data. 

Based  on  a  baseline  database  for  all  controls,  we  developed  a  lesion  load  approach  to  quantifying  those 
abnormal  white  matter  clusters  as  “lesion  load.”  Our  initial  analysis  demonstrated  that  DTI  lesion  load 
is  significantly  correlated  with  patients’  neurocognitive  performance  at  the  acute  stage  (Pearson 
correlation  coefficient  r=-0.87,  p<0.02).  Meanwhile,  we  are  in  the  process  of  analyzing  the  test-retest 
reliability  analysis  for  DTI  baseline  database. 


DTI-VBA  Lesion  Load  vs.  Delayed  Recall 


DTI-VBA  Lesion  Load 


Figure  1.  Correlations  between  DTI  lesion  load  and  patients’  neurocognitive  data.  As  demonstrated  in 
the  figures,  DTI  lesion  load  (both  TBSS  and  VBA  data)  is  correlated  significantly  with  patients’  overall 
SAC  score  and  delayed  recall.  R  squared  values  are  shown  on  each  figure  for  linear  regression. 


e)  Correlation  MR  imaging  data  with  blood  biomarker  data. 


Our  pilot  data  demonstrated  that  blood  biomarker  GFAP  levels  are  significantly  higher  in  patients  with 
intracranial  bleeding  than  those  without  bleeding.  This  finding  has  never  been  reported  before.  It 
implies  that  a  blood  biomarker  might  be  able  to  serve  as  a  screening  tool  for  MRI  inthe  acute  stage. 


Figure  2.  Box-and-whisker  plots  demonstrating  UCH-L1  and  GFAP  concentrations  in  patients  with 
ventricular  hemorrhages  and  hemorrhagic  contusions  and  in  patients  with  non-hemorrhagic  lesions.  The  black 
horizontal  line  in  each  box  represents  the  median,  with  the  boxes  representing  the  interquartile  range. 
Significant  differences  are  indicated  *  (P<  0.05)  or  **  (P  <0.01)  (Mann-Whitney  U-test). 

f)  Resting  state  functional  MRI  (rsfMRI)  analysis  of  mTBI  patients  at  the  acute  stage. 

Our  rsfMRI  analysis  demonstrated  up-regulation  of  patients  functional  status  at  the  acute  stage,  which 
is  consistent  with  the  increased  cerebral  blood  flow  measured  by  ASL  technique. 

Group  ICA  analysis 

The  group  ICA  was  performed  in  two  steps:  1)  The  number  of  associated  voxels  to  the  networks  for  each 
desired  region  and  voxels'  value  were  measured.  Voxels'  value  indicates  voxel  dependency  to  the 
networks.  2)  "Cross-Validation"  was  carried  out  in  order  to  investigate  the  reliability  and  the 
reproducibility  and  to  verify  the  result  of  the  previous  step.  Among  all  investigated  areas  in  DMN  and 
BGN,  only  the  precuneus  and  the  posterior  cingulate  cortex  have  shown  a  difference  between  the  two 
groups  in  both  the  number  and  the  dependency  of  associated  voxels  (see  Figures  3  and  4).  The  cross 
validation  results  support  this  difference  in  the  precuneus  and  the  posterior  cingulate  cortex  (see 
supplement  3  table  2).  The  voxel-wise  analysis  was  performed  on  the  subgroups'  selected  networks.  The 
two-sample  t-test  analysis  revealed  significant  statistical  differences  for  a  p-value  of  0.001  in  the  posterior 
cingulate  cortex  (50  voxels)  and  precuneus  (442  voxels).  Neither  voxel  value  nor  the  number  of  voxels 
reveals  a  difference  in  the  basal  ganglia  network  regions.  To  sum  up,  ICA  indicates  that  there  is  a  group 
difference  between  healthy  subjects  and  patients  in  various  regions.  The  precuneus,  the  posterior 
cingulate  cortex,  and  the  frontal  lobe  show  significant  differences  between  groups  and  these  differences 
are  constant  among  difference  group  ICA  analysis.  This  proves  that  these  regions  could  discriminate 
between  patients  and  healthy  subjects.  In  the  next  step,  individual  analysis  was  performed  to  investigate 
alteration  in  patients’  resting-state  network  compared  to  healthy  control  subjects  at  the  individual  level. 
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Default  Mode  Network  Default  Mode  Network  Basal  Ganglia  Networks 


Figure  3.  Number  of  associated  voxels  to  the  resting  state  networks:  (a)  and  (b)  Defaidt  mode  network  shows 
difference  between  two  groups  in  the  posterior  cingulate  cortex  (a)  and  precuneus  (b),  and  (c)  Basal  ganglia  network 
does  not  show  a  significant  difference. 


Default  Mode  Network  Default  Mode  Network  Basal  Ganglia  Networks 


Figure  4.  The  voxel  dependency  to  resting  state  networks:  (a)  and  (b)  Default  mode  network  shows  a  difference 
between  the  two  groups  in  the  posterior  cingulate  cortex  (a)  and  precuneus  (b),  and  (c).  The  basal  ganglia  network 
shows  none. 
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Individual  analysis 


All  the  statistical  analyses  were  for  group  ICA  were  carried  out  using  all  three  individual  ICA  methods.  A 
subject  in  which  either  BGN  or  DMN  was  not  extracted  properly  was  eliminated  from  further  statistical 
analysis.  Neither  "dual  regression"  nor  '"back-projection",  was  able  to  discriminate  between  the  groups  in 
an  individual  analysis.  However,  the  suggested  method,  spatial  constraint  ICA  using  the  chosen  template 
("A  baseline  for  the  multivariate  comparison  of  resting-state  networks11  ,  Allen,  E.  A.,  and  et  all  (2011), 
demonstrates  an  intriguing  result.  The  two-sampled  t-test  analysis  on  components  extracted  by  using  the 
proposed  method  reveals  a  statistical  significant  difference  between  two  groups  especially  in  the  posterior 
cingulate  cortex  and  surrounding  areas  like  the  precuneus  in  DMN.  For  p-value  =  0.05,  the  statistical  two- 
sample  t-test  map,  which  was  corrected  by  the  spatial  threshold  using  MonteCarlo  simulation,  manifests  a 
cluster  with  137  voxels  which  includes  the  posterior  cingulate  cortex  (55  voxels)  (see  Figure  5). 


Figure  5.  Two-sampled  t-test  results  showing  the  difference  in  default  mode  network  (DMN)  between  two  groups  in 
individual  independent  component  analysis  (ICA).  Panel  (a)  shows  the  cluster  which  is  statistically  significant 
(p=0.05),  and  panel  (b)  shows  the  result  in  the  posterior  cingulate  cortex  (55  voxels).  The  warm  color  labels  the 
region  where  the  healthy  group  has  more  activation  than  patients. 

Seed  point  Analysis 

We  also  performed  seed  point  analysis,  which  is  the  most  common  method  in  functional  connectivity  to 
evaluate  the  alteration  in  larger  scope  using  univariate  approach.  The  activation  map  that  was  extracted 
for  each  subject  using  the  correlation  between  all  voxels  in  the  region  and  whole  brain  voxels  did  not 
show  any  significant  difference  between  healthy  and  patient  groups.  However,  the  statistical  analysis, 
which  used  the  average  of  the  time  series  in  the  region  and  measured  the  group  difference,  revealed 
statistically  significant  results.  In  the  posterior  cingulate  cortex  map,  univariate  functional  connectivity 
analysis  shows  that  the  correlation  map  in  patients  group  is  larger  than  healthy  subjects  group,  and  it  is 
more  distributed  in  patients’  brain  (includes  more  region  of  brain).  Patients  showed  a  more  powerful 
connection  between  the  cingulate  posterior  cortex  and  the  frontal  lobe  regions.  For  a  p-value  of  0.01, 
control  subjects  group  includes  no  functional  connectivity  in  the  dorsolateral  prefrontal  cortex  (BA  9)  and 
the  anterior  cingulate  cortex  (BA  32),  and  it  includes  a  few  voxels  in  the  anterior  prefrontal  cortex  (BA 
10),  while  patient  group  is  significantly  correlated  in  these  regions.  One-sample  t-test  in  seed  point 
analysis  also  demonstrated  a  group  difference  in  PCC  and  precuneus,  which  is  consistent  with  ICA  data. 
The  two-sample  t-test  shows  differences  between  two  groups.  Notably  the  dorsolateral  prefrontal  cortex 
(BA  9)  and  adjoining  voxels  in  Brodmann  area  8,  and  the  anterior  cingulate  cortex  (BA  32)  (see  Figure 
6).  For  p-value=0.01,  the  two-sampled  t-test  map  has  two  clusters  with  117  and  223  voxels.  These 
clusters  include  77  and  87  voxels  in  the  dorsolateral  prefrontal  cortex  (BA  9)  respectively,  and  the  biggest 
cluster  contains  44  voxels  in  the  dorsal  anterior  cingulate  cortex  (BA32)  which  shows  the  importance  of 
these  regions’  liaison  in  mild-TBI  patients. 
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Figure  6.  Two-sample  t-test  (p=0.01)  for  the  posterior  cingulate  cortex  correlation  map:  The  cold  color  labels  the 
region  that  the  patients  group  has  more  correlation  with  the  posterior  cingulate  and  the  warm  color  shows  the  voxels 
that  healthy  group  has  more  correlation  with  the  posterior  cingulate  cortex.  Cross  bar  located  in  different  positions  in 
images  a  and  b. 

In  the  precuneus  map,  the  between  group  comparison  for  p-value  =  0.05  reveals  two  regional  differences. 
The  first  statistically  significant  region  is  in  the  supramarginal  gyrus  with  82  voxels;  the  other  cluster  has 
85  voxels,  which  includes  Brodmann  area  8  and  the  interface  of  the  mid-cingulate  cortex  ,  the  frontal 
superior,  and  the  sup-motor  area  (interface  Brodmann  area  8,  32)  (see  Figure  7.).  Interestingly,  almost  the 
same  area  is  statistically  significant  when  the  posterior  cingulate  cortex  or  the  thalamus  is  a  seed  region. 
For  one-sample  t-test  and  more  details  see  supplements. 


Figure  7.  Two-sample  t-test  (p — 0.05):  The  precuneus  correlation  map  shows  group  difference  in  the  supramarginal 
gyrus  and  the  interface  Brodmann  areas  8  and  32.  The  cold  color  labels  the  region  that  Patients  group  has  more 
correlation  with  the  precuneus.  Additionally,  the  activated  area  in  BA  32  and  BA  8  are  also  activated,  with  spatial 
variation,  when  the  posterior  cingulate  or  the  thalamus  are  used  as  seed  point. 


In  the  next  step,  the  thalamus  correlation  map  was  investigated.  For  p-value  =  0.01,  between  group 
comparison  shows  the  statistically  significant  difference  between  groups  in  the  anterior  prefrontal  cortex 
(BA  10)  (see  Figure  8(a))  and  part  of  the  supramarginal  gyrus  (BA  40)  (see  Figure  8(b)),  which  indicates 
the  effect  of  brain  injury  on  the  thalamus  connections. 
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Figure  8.  Two-sample  t-test  (p=0.01):  The  region  of  interest  is  in  the  thalamus  and  cold  color  shows  the  regions  that 
patient  group  has  more  correlation  with  seed  point  compare  to  healthy  group,  (a)  indicate  that  there  is  a  statistical 
difference  between  two  group  at  the  anterior  prefrontal  cortex  (BA  10)  and  (b)  manifest  the  difference  in  the 
supramarginal  gyrus  (BA  40). 

For  the  amygdala  correlation  map,  the  healthy  group  shows  higher  connectivity  inside  the  amygdala.  For 
the  one-sample  t-test  and  p-value  =  0.01,  healthy  group  includes  73  and  67  voxel  while  for  patients 
group  this  number  decrease  to  39  and  22  (see  Figure  9).  The  two-sample  t-test  (p-value  =  0.01) 
demonstrate  significantly  increased  connectivity  in  the  left  parietal  superior  cortex  in  the  patients  group 
(cluster  size  =  104)  (see  Figure  9). 


Figure  9.  Two-sample  t-test  (p — 0.0 1 )  for  the  amygdala  correlation  map:  The  cold  color  labels  the  region  (the  left 
parietal  superior  cortex  )  that  patients  group  has  more  correlation  with  the  amygdala. 

The  univariate  functional  connectivity  analysis  for  hippocampus  as  ROI  demonstrates  significant 
alteration  in  mTBI  patients.  For  p-value  =  0.01,  three  clusters  (cluster  sizes=  135  ,61,  and  52)  are 
significantly  difference  between  two  groups.  While  the  increased  connectivity  in  the  fusiform  gyrus  and 
the  somatosensory  association  cortex  (BA  7)  are  observed  in  the  patients  group,  the  decreased  intrinsic 
connectivity  is  observed  in  the  BA  48  (see  Figure  10). 


Figure  10.  Two-sample  t-test  (p=0.01)  for  the  hippocampus  correlation  map:  The  cold  color  labels  the  region  that 
the  patients  group  has  more  correlation  with  the  hippocampus  while  hot  color  labels  the  region  that  the  healthy 
group  has  more  correlation.  In  (a)  cross  bar  located  on  BA,  in  (b)  cross  bar  located  on  BA  48,  and  in  (b)  cross  bar 
located  on  the  fusiform  gyrus. 

We  could  not  find  any  statistical  significant  correlation  between  the  neurocognitive  tests,  delayed  recall 
and  SAC  score,  and  any  of  statistical  analysis  results.  The  largest  correlation  value  (0.42)  was  found  in 
the  seed  point  analysis  for  PCC  correlation  map  between  the  mean  of  the  correlation  value  of  the  dorsal 
anterior  cingulate  and  delay  recall  test. 


g)  Develop  susceptibility  weighted  imaging  and  mapping  (SWIM)  results  to  quantify  the  amount  of  iron 
in  microbleeds,  to  monitor  any  changes  (evolution)  of  the  microbleed  over  time  and  to  monitor  oxygen 
saturation  compared  to  normal  volunteers. 

We  have  had  a  very  successful  two  years  in  enhancing  the  ability  of  SWIM  to  quantify  iron  in  the  brain. 
Our  work  has  appeared  in  numerous  papers  listed  at  the  end  of  this  report  and  one  of  our  data  sets  for  both 
SWI  and  SWIM  was  used  as  an  image  on  the  front  cover  of  the  book:  Cerebral  Microbleeds: 
Pathophysiology  to  Clinical  Practice,  edited  by  David  Werring,  published  2012. 

The  main  efforts  we  have  had  in  this  direction  has  been  to  understand  how  to  best  reproduce  quantitative 
susceptibility  maps  given  the  ill-posedness  of  the  solution  using  a  single  orientation  acquisition  approach. 
To  address  the  questions  of  microbleeds,  we  published  a  paper  on  “ Quantitative  susceptibility  mapping  of 
small  objects  using  volume  constraints  ”  (see  ref.  1  below)  where  we  showed  that  even  though  absolute 
susceptibility  may  be  inaccessible  for  small  bleeds,  the  magnetic  moment  could  in  fact  be  well  calculated. 


Bfb 


We  followed  this  work  up  with  a  new  iterative  approach  to  QSM  or  SWIM  as  we  call  it  that  helped  to 
correct  for  some  of  the  current  weaknesses  of  the  methods.  This  paper  was  entitled:  “ Improving 
susceptibility  mapping  using  a  threshold-based  k-space/image  domain  iterative  reconstruction  approach ” 
(see  ref.  2  below).  This  method  was  designed  specifically  to  help  detect  changes  in  oxygen  saturation  that 
might  in  the  future  be  used  to  estimate  perfusion  deficit  in  mTBI.  Although  we  have  not  done  this  study 
yet,  we  believe  that  this  is  a  very  important  direction  for  future  research  and  clinical  care  in  mTBI.  We 
also  published  a  chapter  related  to  mTBI,  imaging  veins,  SWIM  and  perfusion  entitled:  “ The  Presence  of 
Venous  Damage  and  Microbleeds  in  Traumatic  Brain  Injun >  and  the  Potential  Future  Role  of 
Angiographic  and  Perfusion  Magnetic  Resonance  Imaging”  (see  ref.  7  below). 

Still,  there  are  practical  issues  in  getting  the  right  phase  information  and  getting  it  unwrapped  rapidly.  To 
address  this  issue  we  studied  a  multi-echo  approach  to  doing  phase  unwrapping  and  doing  it  on  a  pixel- 
by-pixel  basis  so  that  it  would  be  rapid  and  robust.  This  led  to  another  paper:  “ Catalytic  multiecho  phase 
unwrapping  scheme  (CAMPUS)  in  multiecho  gradient  echo  imaging:  Removing  phase  wraps  on  a  voxel- 
bv-voxel  basis,  ’’(see  ref.  6  below).  Although  this  technique  used  1 1  echoes,  we  are  now  trying  to 
accomplish  this  with  two  echoes  instead. 

Once  we  have  a  good  working  SWIM  dataset,  it  becomes  possible  to  do  something  that  leads  to  a  new 
type  of  SWI  dataset,  what  we  refer  to  as  “True  SWI”  that  uses  the  susceptibility  map  to  create  a  mask  that 
does  not  suffer  from  orientation  effects  on  the  phase  and  is  a  local  mask  again  not  suffering  from  non¬ 
local  field  or  phase  effects.  This  has  been  tentatively  accepted  for  publication  and  is  entitled:  “ True 
susceptibility  weighted  imaging.  Tentatively  accepted  in  J  Magn  Reson  Imaging.  ”  (see  ref.  1  in  the  under 
preparation  section). 

Imaging  iron  and  even  mapping  iron  is  now  possible  with  MRI.  Our  approach  is  straight  forward  and  our 
reconstruction  is  rapid  thanks  to  the  iterative  method  we  propose.  However,  knowing  exactly  how  SWIM 
data  corresponds  with  actual  iron  is  more  difficult.  Along  these  lines,  we  did  a  comparison  of  MRI  with 
XRF  and  published  the  following  paper:  “Measuring  iron  in  the  brain  using  quantitative  susceptibility 
mapping  and  X-ray  fluorescence  imaging.  ”  (see  ref.  3  below). 

Imaging  veins  and  microbleeds  is  important  but  there  is  more  to  the  story  than  just  veins  in  the 
vasculature.  For  more  severe  impacts,  the  arterial  system  may  well  be  compromised  as  well.  In  that  case, 
we  want  to  be  able  to  image  the  vasculature  as  well  as  possible,  and  that  can  now  be  done  with  or  without 
a  contrast  agent  using  our  double  echo  SWI  sequence  as  shown  in:  “ Noncontrast-enhanced  magnetic 
resonance  angiography  and  venography  imaging  with  enhanced  angiography .” 

Finally,  bring  together  other  aspects  of  imaging  along  with  detection  of  hemorrhage  and  the  need  for 
neuroradiologists  to  be  involved,  we  looked  at  combining  DTI  with  SWI  and  other  conventional  imaging. 
Here  we  published  the  paper:  “ Detection  of  hemorrhagic  and  axonal  pathology >  in  mild  traumatic  brain 
injury  using  advanced  MRI:  implications  for  neurorehabilitation." 


In  summary,  we  have  created  an  entirely  new  foundation  for  studying  mTBI  using  SWI  and  SWIM  and 
in  combining  these  with  DTI  and  hopefully  in  the  future  with  perfusion  imaging.  It  is  our  hope  that  this 
will  lead  to  not  only  better  diagnosis  but  hope  for  better  treatment  for  patients. 


KEY  RESEARCH  ACCOMPLISHMENTS 


The  following  accomplishments  have  been  made  in  the  past  year: 

•  Recruitment  and  follow  up  of  patients  and  controls.  Now  we  are  in  the  fast  lane  of  scanning 
subjects  on  average  4+  scans  each  week.  With  the  approval  of  a  no-cost  extension,  we  are  very 
confident  of  recruiting  total  50  patients  and  50  controls  as  our  original  proposal. 

•  Analysis  of  preliminary  data.  Our  data  demonstrated  that,  in  addition  to  an  increased  blood  flow 
as  an  early  response  of  mild  TBI  patients,  patients’  resting  functional  connectivity  tends  to 
increase  as  well. 

•  Blood  biomarker  GFAP  levels  significantly  increased  in  intracranial  cases.  This  implies  a 
screening  tool  of  biomarker  might  be  able  to  serve  for  MRI  identification. 

•  DTI  abnormalities  account  for  patients’  neurocognitive  symptoms  at  the  acute  stage. 

REPORTABLE  OUTCOMES 

Papers  in  preparation  based  on  the  material  discussed  above.  Note  that  in  this  section  and  the  next  section 
that  an  asterisk  means  that  the  paper  has  referenced  the  support  of  this  award  as  will  future  papers. 

a)  *S.  Liu,  K.  Mok,  J  Neelavalli,  YC  Cheng,  J.  Tang,  Y.  Ye  and  E.M.  Haacke,  True 
susceptibility  weighted  imaging.  Tentatively  accepted  in  J  Magn  Reson  Imaging. 

b)  *Zhifeng  Kou,  Ramtilak  Gattu,  Firas  Kobeissy,  Robert  Welch,  Brian  O’Neil,  John  Woodard, 
Sayed  Imran  Ayaz,  Andrew  Kulek,  Robert  Kas’shamoun,  Valerie  Mika,  Conor  Zuk,  E  Mark 
Haacke,  Ronald  Hayes,  Francesco  Tomasello,  Stefania  Mondello.  Combining  Biochemical 
and  Imaging  Markers  to  Improve  Diagnosis  and  Characterization  of  Mild  Traumatic  Brain 
Injury  in  the  Acute  Setting:  Results  from  a  Pilot  Study.  Submitted  to  PLoS  ONE,  in  review. 

c)  Z.  Kou  et  al,  another  paper  on  resting  state  fMRI  is  also  in  preparation  for  submission. 
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Our  group  has  given  a  number  of  presentations  in  national  and  international  meetings.  Over  and  above 
those  already  quoted  for  2011  and  2012  in  the  previous  annual  report,  these  include: 

2013  presentations 

Jie  Yang,  Zhifeng  Kou,  Robert  Dean  Welch,  Randall  Benson,  Ramtilak  Gattu,  Valerie  Mika,  and  E.  Mark 
Haacke.  Neuroimaging  biomarkers  of  mild  traumatic  brain  injury  (mTBI)  and  its  recovery:  A 
preliminary  study  in  acute  setting  30pm.  21st  Annual  ISMRM  Salt  Lake  City,  2013.  Electronic  poster 
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E.  Mark  Haacke  and  Zhifeng  Kou.  “Development  of  MRI  Biomarkers  for  Improved  Diagnosis  of  TBI” 
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Grants,  Honors  and  Awards 

A  major  award  (Seed  Grant  Award)  was  received  by  Dr.  Zhifeng  Kou  from  the  International  Society  for 
Magnetic  Resonance  in  Medicine  (ISMRM)  to  investigate  the  relationship  between  MR  imaging  and 
blood  biomarkers  in  acute  detection  and  outcome  prediction  of  mTBI. 


CONCLUSION 

During  the  last  year  we  have  made  major  progress  in  collecting  the  data  on  both  mTBI  patients  and 
normal  controls.  We  have  demonstrated  already  in  this  group  the  value  of  SWI  and  DTI  and  potential 
correlations  with  neuro-psychological  testing.  This  extension  year  will  allow  us  to  both  complete 
collection  of  the  data  and  answer  a  number  of  key  questions  in  the  quest  for  better  diagnosis  and 
understanding  of  the  pathophysiology  in  mTBI  patients. 
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Abstract: 

Background:  Mild  traumatic  brain  injury  (mTBI)  is  a  significant  public  healthcare  burden 
and  its  diagnosis  remains  a  challenge  in  emergency  department.  Serum  blood 
biomarker-based  tests  and  advanced  magnetic  resonance  imaging  (MRI)  techniques 
have  already  demonstrated  their  potentials  to  improve  brain  injury  detection  even  in 
patients  with  negative  computed  tomography  (CT)  finding.  The  objective  of  this  study 
was  to  determine  the  clinical  value  of  a  combinational  use  of  both  blood  biomarkers 
and  MRI  in  mTBI  detection  and  characterization  in  the  acute  setting  (within  24  hours 
after  injury). 

Methods:  Nine  patients  with  mTBI  were  prospectively  recruited  from  the  emergency 
department.  Serum  samples  were  collected  at  the  time  of  hospital  admission  and  every 
6  hours  up  to  24  hours  post  injury.  Neuronal  (Ubiquitin  C-terminal  Hydrolase-LI  [UCH- 
Llj)  and  glial  (glial  fibrillary  acidic  protein  [GFAP])  biomarker  levels  were  analyzed. 
Advanced  MRI  data  were  acquired  at  9±6.91  hours  after  injury.  Patients' 
neurocognitive  status  was  assessed  by  using  the  Standard  Assessment  of  Concussion 
(SAC)  instrument. 

Results:  The  median  serum  levels  of  UCH-L1  and  GFAP  on  admission  were  increased 
4.9  folds  and  10.6  folds,  respectively,  compared  to  reference  values.  Three  patients 
were  found  with  intracranial  hemorrhages  on  SWI,  and  they  all  had  very  high  GFAP 
levels.  Total  volume  of  brain  white  matter  (WM)  with  abnormal  diffusion  tensor  imaging 
(DTI)  fractional  anisotropy  (FA)  measures  were  negatively  correlated  with  patients' 

SAC  scores,  including  memory.  Both  increased  and  decreased  DTI-FA  values  were 
observed  in  the  same  subjects.  Serum  biomarker  level  was  not  correlated  with  patients' 
DTI  data  nor  SAC  score. 

Conclusions:  Blood  biomarkers  and  advanced  MRI  may  correlate  or  complement  each 
other  in  different  aspects  in  mTBI  detection  and  characterization.  GFAP  might  have 
potential  to  serve  as  a  clinical  screening  tool  for  intracranial  bleeding.  UCH-L1 
complements  MRI  in  injury  detection.  The  patients'  neurocognitive  symptoms  were 
accounted  by  impairment  at  WM  tracts. 
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1  Abstract 

2  Background:  Mild  traumatic  brain  injury  (mTBI)  is  a  significant  public  healthcare  burden 

3  and  its  diagnosis  remains  a  challenge  in  emergency  department.  Serum  blood 

4  biomarker-based  tests  and  advanced  magnetic  resonance  imaging  (MRI)  techniques 

5  have  already  demonstrated  their  potentials  to  improve  brain  injury  detection  even  in 

6  patients  with  negative  computed  tomography  (CT)  finding.  The  objective  of  this  study 

7  was  to  determine  the  clinical  value  of  a  combinational  use  of  both  blood  biomarkers  and 

8  MRI  in  mTBI  detection  and  characterization  in  the  acute  setting  (within  24  hours  after 

9  injury). 

10  Methods:  Nine  patients  with  mTBI  were  prospectively  recruited  from  the  emergency 

11  department.  Serum  samples  were  collected  at  the  time  of  hospital  admission  and  every 

12  6  hours  up  to  24  hours  post  injury.  Neuronal  (Ubiquitin  C-terminal  Hydrolase-LI  [UCH- 

13  LI])  and  glial  (glial  fibrillary  acidic  protein  [GFAP])  biomarker  levels  were  analyzed. 

14  Advanced  MRI  data  were  acquired  at  9±6.91  hours  after  injury.  Patients’  neurocognitive 

15  status  was  assessed  by  using  the  Standard  Assessment  of  Concussion  (SAC) 

16  instrument. 

17  Results:  The  median  serum  levels  of  UCH-L1  and  GFAP  on  admission  were  increased 

18  4.9  folds  and  10.6  folds,  respectively,  compared  to  reference  values.  Three  patients 

19  were  found  with  intracranial  hemorrhages  on  SWI,  and  they  all  had  very  high  GFAP 

20  levels.  Total  volume  of  brain  white  matter  (WM)  with  abnormal  diffusion  tensor  imaging 

21  (DTI)  fractional  anisotropy  (FA)  measures  were  negatively  correlated  with  patients’  SAC 

22  scores,  including  memory.  Both  increased  and  decreased  DTI-FA  values  were  observed 

23  in  the  same  subjects.  Serum  biomarker  level  was  not  correlated  with  patients’  DTI  data 

24  nor  SAC  score. 

25  Conclusions:  Blood  biomarkers  and  advanced  MRI  may  correlate  or  complement  each 

26  other  in  different  aspects  in  mTBI  detection  and  characterization.  GFAP  might  have 
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potential  to  serve  as  a  clinical  screening  tool  tor  intracranial  bleeding.  UCH-L1 
complements  MRI  in  injury  detection.  The  patients’  neurocognitive  symptoms  were 
accounted  by  impairment  at  WM  tracts. 
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1  INTRODUCTION 

2  Traumatic  brain  injury  (TBI)  is  a  significant  public  healthcare  burden  in  the  United  States. 

3  It  accounts  for  1 .7  million  incidents  in  the  United  States  each  year  [1 , 2],  Most  TBI 

4  patients  belong  to  mild  TBI  (mTBI)  severity  due  to  the  improvement  of  motor  vehicle 

5  safety  design  in  recent  years  [3].  mTBI  has  over  1  million  emergency  visits  in  the  United 

6  States  each  year  [4].  It  causes  a  constellation  of  physical,  cognitive,  and  emotional 

7  symptoms  that  significantly  impact  the  patients’  quality  of  life,  and  costs  the  nation  $1 6.7 

8  billion  each  year  [3,  5,  6].  To  date,  most  mTBI  patients  stay  in  the  emergency 

9  department  (ED)  for  only  a  few  hours  and  then  are  discharged  home  without  a  concrete 

10  follow-up  plan.  Up  to  50%  of  mTBI  patients  may  develop  prolonged  neurocognitive 

11  problems  within  one  month  [7,  8],  and  5-1 5%  of  them  continue  to  manifest 

12  neurocognitive  sequelae  at  one  year  [7,  9],  Often,  their  neurocognitive  outcomes 

13  inconsistently  correlate  with  clinical  measures  such  as  the  Glasgow  Coma  Scale  (GCS) 

14  score  and  post-traumatic  amnesia.  Most  mTBI  patients  do  not  present  abnormalities  on 

15  computed  tomography  (CT)  at  admission  and  conventional  magnetic  resonance  imaging 

16  (MRI)  [10,  11]  in  the  emergency  setting.  The  immediate  challenge  for  emergency 

17  physicians  is  to  identify  those  CT-negative  patients  for  intracranial  abnormalities  and 

18  long-term  neurocognitive  symptoms  [12]. 

19 

20  Advanced  MRI  has  demonstrated  improved  sensitivity  in  detecting  TBI  pathologies  and 

21  functional  impairments  that  underlie  patients’  cognitive  symptoms  [13,  14].  Examples 

22  include  diffusion  tensor  imaging  (DTI)  of  axonal  injury  [1 5-24],  susceptibility  weighted 

23  imaging  (SWI)  of  hemorrhagic  lesions  [25-27],  and  others.  Advanced  MRI  offers 

24  anatomical  and  pathological  information,  reflecting  brain  damage  with  high  spatial 

25  resolution  [13,  14].  The  American  College  of  Emergency  Physicians  and  Center  for 

26  Disease  Control  Joint  Study  Panel  highly  recommended  examining  the  role  of  advanced 
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1  MRI  in  the  acute  setting  (within  24  hours  after  injury)  [12].  However,  to  date,  very  few 

2  studies  have  been  designed  for  scanning  mTBI  patients  within  24  hours  after  injury  when 

3  the  patients  are  still  in  ED.  This  is  mainly  due  to  the  high  cost  of  MRI  that  makes  it 

4  prohibitive  to  scan  all  mTBI  patients  and  the  inaccessibility  of  MRI  for  24/7  availability  in 

5  the  acute  setting.  Consequently,  in  the  majority  of  medical  centers  in  North  America, 

6  MRI  is  not  a  standard  of  care  of  mTBI  patients  in  the  acute  setting,  partially  due  to  the 

7  difficulty  of  reimbursement  associated  with  the  high  cost  of  a  MRI  scan. 

8 

9  In  addition  to  patient’s  clinical  characteristics  and  advanced  neuroimaging  studies, 

10  brain-specific  proteins  released  into  biofluids  after  brain  injury,  as  a  result  of  cellular 

11  damage  and  activation,  have  demonstrated  the  potential  to  serve  as  diagnostic  and 

12  prognostic  markers  in  mTBI  [28,  29],  Along  with  providing  improved  diagnostic  capability 

13  and  molecular  characterization  of  subjects  who  sustained  mTBI,  appropriate  biomarker 

14  screening  may  lead  to  a  more  selective  strategy  for  neuroimaging,  reducing  the  need  for 

15  a  substantial  number  of  unnecessary  imaging  exams.  With  these  aims  in  mind,  recently, 

16  several  groups  reported  the  application  of  neuroproteomics  to  identify  and  characterize 

17  biochemical  markers  of  TBI  [30].  Among  the  variety  of  blood  proteins  that  have  been 

18  used  to  investigate  the  neuronal  marker,  Ubiquitin  C-terminal  Hydrolase-LI  (UCH-L1) 

19  [31-33]  and  the  astrocyte-specific  protein  Glial  Fibrillary  Acidic  Protein  (GFAP)  [34],  and 

20  all-Spectrum  Breakdown  Products  (SBDPs)  for  axonal  injury  [35,  36]  seem  particularly 

21  promising. 

22 

23  UCH-L1  is  a  small  (25  kDa),  neuronal  protease  involved  in  either  the  addition  or  removal 

24  of  ubiquitin  from  proteins  that  are  destined  for  metabolism  via  the  ATP-dependent 

25  proteasome  pathway;  it  is  highly  enriched  in  the  brain  (1  -  5%  of  total  soluble  brain 

26  protein)[37,  38].  Mutations  and  polymorphisms  of  UCH-L1  have  been  associated  with 
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1  familial  Parkinson’s  Disease  [39].  UCH-L1  is  released  into  the  extracellular  space  as  a 

2  consequence  of  cell  destruction  under  diverse  pathological  conditions  affecting  the  brain. 

3  Previous  clinical  studies  have  demonstrated  increased  UCH-L1  levels  in  cerebral  spinal 

4  fluid  (CSF)  and  in  serum  in  severe  TBI  patients  and  that  the  magnitude  of  this  increase 

5  correlated  with  injury  severity,  CT  finding  and  patient  outcome  [31 , 33].  Recently,  a  study 

6  was  completed  investigating  UCH-L1  in  adults  with  mild  and  moderate  TBI  showing 

7  increased  UCH-L1  levels  in  mTBI  patients  compared  to  uninjured  controls  and  that  UCH- 

8  LI  was  able  to  detect  intracranial  lesions  on  CT  with  an  area  under  the  curve  (AUC)  of 

9  0.73  [40].  Based  on  these  encouraging  results  and  the  fact  that  UCH-L1  is  specific  to 

10  neurons  and  its  high  specificity  and  abundance  in  the  CNS,  it  appears  to  be  an  excellent 

11  candidate  biomarker  for  the  brain  injury  clinical  studies. 

12 

13  Compared  with  UCH-L1 ,  GFAP  is  a  monomeric  intermediate  filament  protein,  major 

14  constituent  of  the  astroglial  cytoskeleton,  and  highly  brain-specific  [41 , 42],  This  glial 

15  protein  represents  an  ideal  complement  of  the  neuronal  UCFI-L1 ,  as  demonstrated  by  a 

16  recent  study  showing  that  the  correlations  between  these  2  markers  reflect  structural 

17  changes  detected  by  neuroimaging  and  may  be  used  as  an  indicator  for  differing 

18  intracranial  pathologies  after  brain  trauma  [43].  Additionally,  previous  studies  evaluating 

19  severe  TBI  patients  demonstrated  that  GFAP  concentrations  were  associated  with  injury 

20  severity  and  outcome  [44,  45].  Recently,  in  a  prospective  cohort  study  of  108  patients 

21  with  mild  or  moderate  TBI,  GFAP  was  found  to  be  elevated  in  the  serum  within  1  h  after 

22  injury,  discerning  TBI  patients  from  uninjured  controls  with  an  area  under  the  curve 

23  (AUC)  of  0.90  and  discriminating  patients  with  and  without  intracranial  lesion  as 

24  assessed  by  CT  with  an  AUC  of  0.79  [46]. 

25 
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1  As  a  simple  biofluid-based  rapid  diagnostic  tool,  serum  biochemical  markers  offer  great 

2  potential  for  rapid,  accurate,  and  cost-effective  diagnosis  of  brain  injury,  and  a  temporal 

3  profile  of  blood  protein  levels  might  be  indicative  of  disease  progression  or  resolution. 

4  On  the  other  hand,  emerging  data  suggest  that  MRI  and  especially  advanced  MRI 

5  techniques  are  very  sensitive  in  detecting  brain  injury  that  are  occult  in  clinical  imaging 

6  by  providing  spatial  and  pathophysiological  information.  Although  the  combinational  and 

7  complementary  use  of  these  tools  is  promising  and  might  have  important  implications  for 

8  improving  injury  detection  and  outcome  prediction,  the  correlation  among  injury 

9  pathologies  at  tissue  level  assessed  by  neuroimaging  and  at  the  protein  level  as 

10  assessed  by  blood  biomarker  profiles  has  not  yet  been  elucidated. 

11 

12  Our  objective  in  this  study  was  to  evaluate  serum  GFAP  and  UCH-L1  levels  after  mTBI 

13  and  their  correlation  to  the  advanced  MRI  findings  in  a  pilot  cohort  in  an  acute  setting.  In 

14  particular,  we  were  interested  in  establishing  a  serum  profile  of  these  biomarkers  that 

15  might  serve  as  signatures  for  the  presence  of  brain  pathology  as  assessed  by  advanced 

16  MRI  methods,  and  thus  aid  in  the  identification  of  patients  who  need  an  MRI  scan  in  the 

17  acute  setting. 

18 

19  MATERIALS  AND  METHODS 

20  Patient  Recruitment 

21  This  study  was  approved  by  both  the  Human  Investigation  Committee  of  Wayne  State 

22  University  and  the  Institutional  Review  Board  of  Detroit  Medical  Center.  Written  informed 

23  consent  was  obtained  from  each  subject  before  enrollment. 

24 

25  A  total  of  9  patients  who  sustained  mTBI  were  prospectively  recruited  from  the  ED  of 

26  Detroit  Receiving  Hospital  (DRH),  a  Level-1  trauma  center,  which  is  an  affiliated  hospital 
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1  of  Detroit  Medical  Center  (DMC).  Patient  eligibility  was  based  on  the  mTBI  definition  by 

2  the  American  Congress  of  Rehabilitation  Medicine  [47]  with  the  following  inclusion 

3  criteria:  Patients  aged  18  or  older  with  an  initial  Glasgow  Coma  Scale  (GCS)  score  of 

4  13-15  in  ED  with  any  period  of  loss  of  consciousness  less  than  30  minutes  or  any  post 

5  traumatic  amnesia  less  than  24  hours,  or  recorded  change  of  mental  status  (confused, 

6  disoriented  or  dazed).  All  patients  required  a  CT  scan  as  part  of  their  clinical  evaluation. 

7  All  of  them  were  be  able  to  speak  English.  The  exclusion  criteria  included  patients  under 

8  the  age  of  1 8  years,  pregnant  woman,  medically  documented  history  of  brain  injury, 

9  neurological  disorders  or  psychoactive  medications,  history  of  substance  abuse,  CT 

10  indication  of  any  metal  in  the  brain  and  body,  known  contraindication  to  MRI  such  as  a 

11  pacemaker  or  other  non-MR  compatible  implanted  device  as  defined  by  metal  screening 

12  procedure,  or  patients  without  a  clear  history  of  trauma  as  their  primary  event  (e.g., 

13  seizure,  epilepsy,  etc).  In  the  acute  stage,  a  patient  might  have  mental  status  change  or 

14  amnesia  in  which  medical  history  may  not  be  properly  obtained,  thus  the  patient's  record 

15  was  retrospectively  screened  as  well  to  exclude  any  patient  who  does  not  fit  our 

16  inclusion  criteria.  Additionally,  we  performed  an  imaging  study  of  18  healthy  controls 

17  without  history  of  head  injury  or  antecedents  of  central  nervous  system  disease. 

18 

19  Neurocognitive  Assessments 

20  At  the  acute  setting,  once  a  patient  was  conscious  and  stable,  they  were  administered 

21  neurocognitive  tests  and  surveyed  about  their  post-concussion  symptoms  (PCS).  Given 

22  the  situation  of  emergency  care,  it  is  not  feasible  to  perform  a  full  battery  of 

23  neuropsychological  assessment.  Instead,  a  short  instrument  called  Standardized 

24  Assessment  of  Concussion  (SAC)  [48]  was  used  to  assess  the  patients’  neurocognitive 

25  status.  The  SAC  instrument  was  originally  developed  for  onsite  testing  of  subject’s 

26  neurocognitive  performance  after  sports  concussion  [49],  It  has  been  reported  that  SAC 
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1  is  sensitive  to  the  acute  changes  following  concussion  and  it  only  requires  limited 

2  training  of  an  administrator  [50].  The  SAC  assesses  4  cognitive  domains  including 

3  orientation,  attention,  immediate  memory  and  delayed  recall,  and  the  resulting  points 

4  give  a  patient  score  between  0  (indicating  greater  cognitive  deficit)  and  30.  Previous 

5  results  report  its  sensitivity  to  brain  injury  in  the  emergency  setting,  particularly  in  that 

6  delayed  recall  is  more  sensitive  to  brain  injury  [50].  The  Emergency  Room  Edition  of  the 

7  SAC  instrument  also  has  a  graded  symptom  checklist  with  all  PCS  symptoms  listed.  The 

8  patients  were  asked  to  grade  each  symptom  from  none,  mild,  moderate,  to  severe, 

9  (graded  from  a  0  to  3  respectively).  The  total  points  were  the  overall  PCS  score. 

10 

11  Neuroimaging  Protocol 

12  In  the  ED,  once  a  patient  was  cleared  of  any  immediate  life-  threatening  risk  after  a  CT 

13  scan  was  completed  and  was  stable  enough  for  an  MRI,  the  patient  was  transported  into 

14  the  MRI  center  for  imaging  scan.  All  MRI  data  were  collected  on  a  3-Tesla  Siemens 

15  Verio  scanner  with  a  32-channel  radio  frequency  head  coil.  Subjects'  head  was  fixed  by 

16  a  foam  pad  to  restrict  motion.  Imaging  protocol  includes  SWI,  DTI,  and  resting  state 

17  functional  MRI,  in  addition  to  the  baseline  structural  imaging  (T1,  T2  gradient  recalled 

18  echo  [GRE]  and  T2  fluid-attenuated  inversion  recovery  [FLAIR])  sequences,  with  total 

19  data  acquisition  time  of  39  min.  SWI  is  a  3-dimentional,  T2*  based  GRE  sequence  with 

20  long  TE  and  3-D  flow  compensation.  The  phase  images  were  high-pass  filtered  (96x96 

21  filter  size)  by  using  an  in-line  manufacture-applied  filter  and  then  integrated  with 

22  magnitude  images  to  generate  the  processed  SWI  image  to  better  delineate  the  spatial 

23  relation  between  microhemorrhages  and  veins  [13,  51].  SWI  parameters  include  TR/TE 

24  of  30/20ms,  Flip  angle  of  15  degree,  bandwidth  of  100  Flz/Px,  field  of  view  (FOV)  of 

25  256x256  mm2,  imaging  matrix  of  512x256,  25%  oversampling,  slice  thickness  of  2  mm, 
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1  total  64  slices,  20%  distance  factor,  GRAPPA  iPat  factor  of  2,  with  resultant  voxel  size  of 

2  0.5x1  x2  mm3  and  imaging  acquisition  time  4m  and  18s. 

3 

4  DTI  sequence  is  a  standard  echo  planar  imaging  (EPI)  2D  sequence  provided  as  part  of 

5  the  Siemens  package.  The  parameters  include  TR/TE  of  13300/124  ms,  EPI  factor  of 

6  192,  bandwidth  of  1240  Hz/Px,  FOV  of  256x256  mm2,  imaging  matrix  of  192x192,  slice 

7  thickness  of  2  mm,  total  60  slices,  no  gap  between  slices,  30  gradient  directions,  2 

8  averages,  B  values  of  0/1000  s/mm2,  anterior-posterior  phase  encoding,  GRAPPA 

9  acceleration  factor  of  2,  with  resultant  voxel  size  of  1 .3x1 .3x2  mm3  and  data  acquisition 

10  time  of  14m  26s. 

11 

12  Image  Processing  and  Interpretation 

13  All  SWI  and  DTI  images  were  processed  by  a  co-author  who  was  blinded  to  the  subjects’ 

14  clinical  conditions  to  avoid  any  bias.  All  SWI  images  were  further  processed  by  using  our 

15  in-house  software  SPIN  (signal  processing  for  NMRI)  (MRI  Research  Institute,  Detroit, 

16  Michigan).  All  structural  MRI  images,  including  SWI  images,  were  read  by  two  board 

17  certified  neuroradiologists  to  identify  other  conditions  that  may  confound  the  findings. 

18  The  neuroradiogists  were  blinded  to  the  medical  history  and  conditions  of  subjects  to 

19  avoid  any  bias  as  well.  We  also  graded  the  structural  imaging  findings  based  on  their 

20  radiologic  report  for  statistical  analysis:  0  for  negative  finding;  1  for  non-specific  finding, 

21  including  non-specific  WM  hyperintensities;  and  2  for  traumatic  hemorrhage. 

22 

23  DTI  image  processing:  Preprocessing  of  DTI  images  was  carried  out  by  using  DTI 

24  Studio  (https://www.dtistudio.org).  The  preprocessing  steps  included  motion  correction 

25  and  eddy  current  correction  by  using  automatic  image  registration  (AIR)  of  all  the 

26  diffusion  weighted  images  to  B0  image  as  a  reference.  Fractional  Anisotropy  (FA)  maps 
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1  were  created  from  the  tensor  calculations  by  suppressing  the  background  noise  on  BO 

2  image  in  DTI  Studio.  Skull  stripping  of  the  BO  image  was  done  by  using  the  BET  routine 

3  package  in  Mricro  with  a  fractional  intensity  of  0.1 .  A  binary  mask  of  the  skull-stripped  BO 

4  image  was  used  on  the  FA  map  to  preserve  the  pure  brain  parenchyma. 

5 

6  Voxel  based  analysis  (VBA):  VBA  was  performed  after  the  skull  stripped  FA  image  was 

7  spatially  normalized,  by  using  a  non-linear  algorithm,  to  the  standard  FMRIB  FA 

8  template  for  all  the  27  subjects  (9  mTBI  patients  and  18  healthy  controls)  in  SPM8 

9  software 

10  (http://www.fil.ion.ucl.ac.uk/spm/software/spm8/).  A  mean  FA  map  and  a  standard 

11  deviation  of  FA  map  were  created  for  the  controls.  A  Z-score  map  was  created  for  each 

12  individual  patient  in  a  voxel-based  approach:  For  each  voxel,  the  FA  difference  between 

13  each  individual  patient  and  the  mean  of  controls  was  divided  by  the  standard  deviation  of 

14  the  controls  at  the  same  voxel.  Voxels  with  Z-score  >  2  for  increased  FA  and  Z-score  <  - 

15  2  for  decreased  FA  were  selected  after  thresholding  as  abnormal  for  further 

16  consideration. 

17 

18  Tract  based  spatial  statistics  ( TBSS ):  Similar  approach  was  also  used  to  evaluate  the 

19  lesion  load  by  using  TBSS  method  after  registering  each  subject  nonlinearly  to  the 

20  standard  FMRIB  FA  template.  TBSS  analytical  approach  was  used  to  compare  patient 

21  and  control  groups  to  evaluate  the  WM  changes  after  TBI.  All  the  processing  steps  were 

22  performed  according  to  the  TBSS  manual 

23  (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/TBSS/UserGuide).  Briefly,  every  subject’s  FA  image 

24  was  spatially  and  non-linearly  normalized  to  a  standard  FA  FMRIB  template  and 

25  transformed  into  a  standard  space  using  FNIRT  algorithm  from  FSL  software  package. 

26  Subsequently,  a  mean  FA  image  was  created  from  this  set  of  non-linearly  transformed 
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1  images.  A  search  algorithm  then  created  a  mean  skeleton,  looking  for  the  local  maxima 

2  perpendicular  to  the  WM  tract  across  the  whole  brain  volume  in  all  the  transformed 

3  images,  and  then  projected  this  skeleton  across  all  the  subjects  in  the  group  to  extract 

4  skeleton  of  individual  subject.  Then  a  voxel-wise  permutation-inference  analysis  was 

5  carried  out  between  the  skeletons  of  two  groups,  and  a  two  tailed  t-statistics  was 

6  performed  to  extract  the  voxels  that  fall  below  or  above  a  certain  threshold.  These  voxels 

7  were  converted  to  a  p  value  based  on  the  threshold  set  by  the  t-stat  value  and  the 

8  cluster  size.  Parameters  used  for  the  TBSS  analysis  were  a  skeleton  threshold  of  0.3  to 

9  eliminate  grey  matter  (GM)  voxels  or  partial  volume  effects  and  cluster  forming  threshold 

10  t  of  4. 

11 

12  Masking  out  non-white  matter  voxels:  The  selected  abnormal  voxels  were  further  filtered 

13  to  eliminate  the  GM  and  the  non-WM  voxels  accounting  for  the  partial  volume  effects 

14  arising  from  CSF  and  GM.  This  filtering  was  done  by  segmenting  each  non-  linearly 

15  normalized  FA  map  and  the  mean  FA  map  from  the  controls  into  WM,  GM  and  CSF  in 

16  SPM8  and  consequently  creating  a  compound  mask  by  applying  a  threshold  (p>0.78)  on 

17  these  two  segmented  WM  FA  images.  In  this  way,  we  were  able  to  get  rid  of  the  false 

18  positives  arising  on  the  edges  of  WM  because  of  the  mis-registration,  and  voxels  having 

19  FA  <  0.3  were  considered  to  be  non  WM  voxels  and  discarded  by  using  a  mask. 

20  Spurious  voxels  that  doesn’t  form  a  cluster  size  of  at  least  1 0  voxels  were  discarded  as 

21  random  noise  in  VBA  and  cluster  size  of  less  than  5  were  discarded  as  random  noise  in 

22  TBSS.  Clusters  were  extracted  using  cluster  tool  in  FSL.  As  a  result,  the  total  number  of 

23  selected  voxels  was  defined  as  lesion  load  for  each  subject  in  both  VBA  and  TBSS 

24  analyses. 

25 

26  Blood  Collection  and  Biomarker  Analysis 
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1  All  blood  samples  were  collected  within  6  hours  after  injury,  upon  subject’s  arrival  to  ED 

2  and  then  every  6  hours  thereafter,  until  discharge  or  up  to  24  hours.  Samples  were 

3  immediately  centrifuged  at  4000  rpm  for  1 0  min  and  frozen  and  stored  at  -80 °C  until  the 

4  time  of  analysis.  Blinded  sample  analysis  was  conducted  in  a  central  laboratory 

5  (Banyan  Biomarkers,  Alachua,  FL)  employing  electro-chemiluminescent  immunoassay 

6  method  (ECL-IA)  for  quantitative  analysis  of  UCH-L1  and  GFAP  in  human  serum 

7  samples  using  the  MSD  platform  (MesoScale  Discovery,  Gaithersburg,  MD).  The  UCH- 

8  LI  assay  system  utilizes  a  mouse  monoclonal  IgM  anti-human  UCFI-L1  antibody  for  solid 

9  phase  immobilization  to  capture  UCF1-L1  from  samples.  The  UCFi-LI  antigen  in  turn 

10  binds  to  a  sulfo-tag  labeled  anti-mouse  antibody.  The  GFAP  ECL-IA  utilizes  a  proprietary 

11  mouse  monoclonal  IgG  anti-human  GFAP  antibody  for  solid  phase  immobilization  and  a 

12  proprietary  polyclonal  rabbit  antibody  for  detection.  The  rabbit  IgG  polyconal  detection 

13  antibody  in  turn  binds  to  a  sulfo-tag  labeled  anti-rabbit  antibody.  Detection  signal  occurs 

14  when  an  electrical  current  is  applied  to  the  electrodes  at  the  bottom  of  each  well  of  the 

15  plate.  The  signal  is  measured  at  620  nm.  Quantitative  determination  of  the  biomarker 

16  concentration  is  achieved  by  comparing  the  unknown  sample  signal  intensities  to  a 

17  standard  curve,  obtained  from  the  calibrators  run  in  the  same  assay.  Target 

18  concentrations  are  reported  in  ng/ml.  Each  assay  plate  included  3  QC  controls  at  high, 

19  medium  and  low  concentrations  of  the  assay  range,  each  plated  in  duplicate.  Calibrators 

20  were  prepared  in  Pooled  Pluman  Serum  (PHS)  matrix.  Specifically,  a  serial  dilution  of  the 

21  calibrator  protein  is  prepared  and  aliquots  of  that  calibrator  solution  are  assayed  in  the 

22  same  assay  volume  and  under  the  same  conditions  as  the  samples.  The  calibrator 

23  signal  intensities  were  used  to  generate  a  dose  response  curve  and  to  calculate  the 

24  sample  concentrations  using  a  weighted  four  parameter  logistic  function  (MSD  software 

25  and  MSD  reader).  The  lower  limit  of  detection  of  the  UCFI-LI  and  GFAP  assays  was 

26  determined  to  be  0.10  and  0.008  ng/mL,  respectively.  Samples  with  undetectable  levels 
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1  of  UCH-L1  or  GFAP  were  assigned  a  value  of  50%  of  the  lower  limit  of  detection  (ie, 

2  0.05  and  0.004  ng/mL,  respectively).  The  median  (IQR)  serum  UCH-L1  and  GFAP 

3  concentrations  determined  in  blood  samples  from  29  healthy  volunteers  using  these 

4  assays  were  used  as  normal  reference  values  (Table  2). 

5 

6  Statistical  analysis 

7  Statistical  analyses  were  performed  by  using  SAS  version  9.2  (Cary,  NC,  USA)and  R 

8  software  (http://www.r-project.org).  Data  normality  was  assessed  by  using  the 

9  Kolmogorov-Smirnov  test.  Results  for  continuous  variables  are  presented  as  mean  (SD) 

10  or  median  (interquartile  range)  as  appropriate.  Frequencies  and  percentages  are 

11  presented  for  categorical  variables.  Between-group  differences  were  assessed  by  the 

12  Student’s  t-test  (for  normally  distributed  continuous  variables)  and  the  Mann-Whitney  U 

13  test  (non-normal  continuous  variables),  Pearson's  chi-squared  test  was  used  to  explore 

14  the  relationships  between  categorical  variables.  Pearson  correlations  were  performed  to 

15  determine  the  relationships  among  different  parameters,  including  imaging,  biomarkers 

16  and  patients  neurocognitive  measurements.  The  relationship  between  biomarker 

17  concentration  and  parameters  for  TBI  severity,  neuroimaging  and  neurocognitive  scores 

18  was  assessed  by  bivariate  correlations  (Spearman’s).  Two  sided  tests  were  used  and 

19  p<0.05  was  considered  significant. 

20 

21  RESULTS 

22  Characteristics  of  the  Subjects 

23  Individual  patient  demographic  and  clinical  characteristics  are  presented  in  Table  1 .  In 

24  our  mTBI  cohort  (n  =  9),  8  (89%)  subjects  were  males  and  1  (11%)  female,  the  average 

25  patient  age  was  41 ,22±14.37  (mean  ±  standard  deviation)  years,  and  all  their  GCS  score 

26  were  15  upon  ER  entrance.  Five  (55%)  patients  were  injured  in  assault  and  4  (45%) 
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1  were  victims  of  motor  vehicle  accidents.  The  median  SAC  and  PCS  scores  were  23.5 

2  and  15.5,  respectively.  Two  patients  presented  positive  findings  in  acute  CT  scan:  one 

3  with  epidural  hematoma  and  cortical  contusion  in  parieto-temporal  region  and  the  other 

4  with  small  subarachnoid  hemorrhage  (SAH).  MRI  scan  were  performed  at  9  ±6.91  hours 

5  after  injury.  Three  patients,  including  those  two  CT  positive  patients,  presented 

6  hemorrhagic  findings  on  structural  MRI.  One  of  these  patients  (Case  1)  presented  small 

7  hemorrhages  that  were  completely  missed  by  CT,  and  another  patient  was  CT  positive 

8  but  missed  the  intraventricular  hemorrhage  by  CT.  In  the  control  population  (n=18),  61% 

9  were  men  and  39%  women,  and  the  average  patient  age  was  34.83±14.30  years.  There 

10  was  no  age  difference  between  patient  and  controls,  but  a  significant  gender  difference 

11  was  found  between  these  two  groups  (p=0.002,  Chi-square  test). 

12 

13  Patients’  Neurocognitive  Performance 

14  The  mean  patients’  SAC  score  was  22.75  ±SD  2.6.  We  compared  this  mean  with 

15  published  normative  data  of  over  568  subjects  [52]  (mean  ±  SD,  26.3±2.2).  The 

16  patients’  mean  SAC  score  was  significantly  below  this  published  mean  score  (t(8)=- 

17  4.180,  p=0.0041).  Among  each  subcategory  of  SAC  test,  patients’  delayed  recall  and 

18  immediate  memory  were  both  significantly  lower  than  published  normalized  data 

19  (p=0.042  and  p=0.021 ,  respectively),  and  concentration  had  a  trend  towards  significance 

20  (p=0.056). 

21 

22  MRI  Findings 

23  Both  increased  and  decreased  FA  beyond  the  threshold  (t>=4)  were  found  in  all  patients 

24  with  variable  number  of  clusters  in  different  locations  of  the  WM.  By  adding  these 

25  clusters  together,  the  volume  of  abnormal  FA,  called  “lesion  load”,  was  used  to  correlate 

26  with  patients’  neurocognitive  and  biomarker  data.  DTI  “lesion  load”  with  pure  increased 
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1  FA  was  found  significantly  higher  than  that  with  pure  decreased  FA  (student  t-test, 

2  p=0.034  for  TBSS  analysis;  and  p=0.017  for  VBA  analysis).  This  suggests  the  increased 

3  FA  as  the  main  pathology. 

4 

5  Figure  1  presents  the  relationship  between  SAC  scores  and  overall  DTI  lesion  load, 

6  which  contains  total  number  of  voxels  with  abnormal  FA  (either  increased  or  decreased). 

7  Specifically,  SAC  scores  were  found  to  be  inversely  correlated  with  DTI-TBSS  lesion 

8  load  (Pearson  r=-0.883,  p  =0.004)  and  DTI-VBA  lesion  load  (r=-0.796,  p=0.018)  and  had 

9  an  almost  significant  correlation  with  age  (r=-0.701 ,  p  =0.053).  For  subcategory  of  SAC 

10  scores,  both  DTI-TBSS  lesion  load  and  DTI-VBA  lesion  load  were  correlated  with  SAC 

11  delayed  recall  (r=-0.834,  p=0.010  and  r=-0.796,  p=0.018,  respectively).  DTI-TBSS  and 

12  DTI-VBA  lesion  loads  were  strongly  correlated  (Pearson  r=0.881  ,p=0.002).  There  was 

13  also  a  partial  correlation  of  SAC  scores  with  DTI-TBSS  and  DTI-VBA  lesion  loads  after 

14  controlling  for  age  (r=-0.893,  p=0.0001 ;  r=-0.82,  p=0.0016,  respectively).There  was  also 

15  a  partial  correlation  of  SAC  delayed  recall  with  patients’  DTI-TBSS  and  VBA  lesion 

16  loads,  controlling  for  age  (r=-0.858,  p=0.001 ;  r=-0.81 1 ,  p=0.002,  respectively). Similarly, 

17  after  controlling  for  age,  DTI-TBSS  and  DTI-VBA  lesion  loads  remained  significantly 

18  correlated  (Pearson  r=0.853  and  p=0.001). 

19 

20  By  further  investigating  the  FA  increase  vs.  FA  decrease  in  association  with  patients’ 

21  SAC  score,  the  lesion  load  with  pure  FA  increase  was  found  significantly  correlated  with 

22  patients  SAC  scores  (Pearson  r=-0.87,  p=0.005  for  TBSS;  and  r=-0.86,  p=0.005  for 

23  VBA)  and  delayed  recall  (Pearson  r=-0.79,  p=0.018  for  TBSS  and  r=-0.770,  p=0.025  for 

24  VBA).  In  contrast,  the  lesion  volume  of  pure  FA  decrease  was  neither  correlated  with 

25  SAC  nor  delayed  recall  (all  p>0.05  for  Pearson  correlation). 

26 
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1  The  structural  imaging  finding,  including  SWI,  was  neither  correlated  with  patients’  SAC 

2  scores  nor  PCS  scores  (all  p>0.05). 

3 

4  Effect  of  non-specific  MRI  findings 

5  Given  the  fact  that  over  half  cases  have  non-specific  findings  on  structural  MRI, 

6  particularly  WM  hyperintensities,  student  T-tests  were  further  performed  to  see  its  effect. 

7  No  group  difference  was  found  between  cases  with  non-specific  findings  and  cases 

8  without  non-specific  findings  (all  with  p>0.3),  in  terms  of  their  SAC  score,  DTI  lesion 

9  loads  (both  TBSS  and  VBA  data),  and  biomarker  levels  (both  UCH-L1  and  GFAP). 

10  Serum  Concentrations  of  UCH-L1  and  GFAP 

11  Median  serum  concentrations  taken  at  the  time  of  hospital  admission  in  the  patients, 

12  within  6  hours  after  injury,  were  raised  4.9  folds  for  UCH-L1 ,  and  1 0.6  folds  for  GFAP 

13  compared  to  the  laboratory  reference  values  in  controls  (Table  2  and  Figure  2).  Serum 

14  UCPI-L1  concentrations  on  admission  did  not  correlate  with  GFAP  (r=-0.24,  p=0.52). 

15  Serum  biomarker  concentrations  at  the  time  of  hospital  admission  did  not  correlate  with 

16  age,  time  to  sample  withdrawal,  GCS,  SAC  and  PCS  score  (p>0.5).  Patients  injured  in 

17  assault  had  significantly  higher  UCFi-LI  concentrations  than  patients  injured  in  a  MVA 

18  (median  0.35  vs  0.10  ng/ml,  p=0.03)  (Figure  3)  while  GFAP  concentrations  were  not 

19  associated  with  mechanism  of  injury. 

20 

21  Neither  UCFI-LI  nor  GFAP  concentrations  were  associated  with  structural  MRI  grading 

22  or  DTI  lesion  load,  as  assessed  by  Pearson  correlation  (p>0.05).  Plowever,  patients  with 

23  hemorrhages  on  structural  MRI  presented  significantly  higher  levels  of  GFAP  compared 

24  with  the  non-  hemorrhagic  group  (p=0.002)  (Table  2  and  Figure  3). 

25 
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1  Temporal  profile  of  biomarker  levels  indicates  that  UCH-L1  tend  to  peak  at  the 

2  admission  (within  6  hours  after  injury)  and  GFAP  at  12  hours  after  injury.  See 

3  supplement  Figure  SI  for  details. 

4  Illustrative  Cases 

5  Case  1  -  Intraventricular  Hemorrhage  Missed  by  CT 

6  A  56-year  old  male  driver  suffered  mental  status  change  after  his  car  was  rear-ended  by 

7  another  vehicle.  Fie  presented  in  the  ED  with  a  GCS  score  of  15  with  abrasion  and  a 

8  small  laceration  on  his  left  eyebrow  without  closure  and  left  clavicle  fracture.  His  major 

9  clinical  symptoms  were  left  shoulder  pain  and  headache.  Non-contrast  CT  scan  showed 

10  no  intracranial  abnormalities.  Initial  MRI  scanning  performed  at  20  hours  after  injury 

11  revealed  small  foci  of  intraventricular  blood  on  the  left  side  and  small  blood  product  in 

12  the  left  lingual  gyrus  and  several  nonspecific  WM  hyperintensies  (Figure  4).  Graphs 

13  displaying  time  course  of  UCF1-L1  and  GFAP  are  shown  in  Figure  4.  In  the  sample 

14  obtained  on  admission,  GFAP  levels  were  markedly  high.  GFAP  elevation  persisted 

15  throughout  the  monitoring  time  gradually  decreasing  at  24  hours  post  injury  (median 

16  4.610,  range  3.241-6.475  ng/ml).  In  contrast,  in  the  same  blood  samples,  UCPI-L1  levels 

17  were  only  slightly  higher  compared  to  controls  (median  0.098,  range  0.055-0.1410 

18  ng/ml). 

19 

20  Case  2  -  Traumatic  axonal  Injury  case  with  normal  structural  MRI  and  high  UCH-L1 

21  level 

22  This  64-year  old  male  patient  was  a  victim  of  an  assault  and  suffered  brief  loss  of 

23  consciousness  and  femur  fracture.  Fie  presented  in  the  ED  with  a  GCS  score  of  1 5  with 

24  symptoms  of  severe  headache,  dizziness,  not  feeling  sharp,  memory  problems,  poor 
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1  concentration,  fatigue/sluggish,  sadness/depression,  and  irritability.  His  SAC  score  was 

2  19  out  of  30,  and  delayed  recall  0  out  of  5.  Non-contrast  CT  scan  showed  no  intracranial 

3  abnormalities.  MRI  scan  performed  at  7  hours  after  injury  demonstrated  multi-foci  of 

4  non-specific  WM  hyper-intensities  on  FLAIR,  but  no  intracranial  bleeding  on  SWI.  DTI 

5  data  revealed  multi-clusters  with  significantly  decreased  FA  in  superior  corona  radiator, 

6  corticospinal  tract  (see  Figure  5),  in  suggestion  of  traumatic  axonal  injury  at 

7  microstructural  level.  The  results  of  neurochemical  of  biomarker  time  course 

8  demonstrated  UCH-L1  being  consistently  elevated  (>0.2ng/ml)  with  high  initial  values 

9  that  gradually  decreased.  GFAP  was  slightly  increased  with  a  peak  at  12  hours  post- 

10  injury  (peak  0.063  ng/ml). 

11 

12  DISCUSSION 

13  To  the  best  of  our  knowledge,  this  is  the  first  effort  to  combine  both  blood  biomarker  and 

14  advanced  MRI  to  improve  the  detection  and  characterization  of  mild  TBI  in  the  acute 

15  setting  (within  24  hours  after  injury).  We  found  that  a)  the  biomarker  levels  were 

16  significantly  higher  in  mTBI  patients  after  injury;  b)  the  levels  of  GFAP  were  highest  in  all 

17  subjects  with  intracranial  bleeding  on  SWI,  which  is  new  finding  in  mTBI  research;  c)  the 

18  total  volume  of  WM  voxels  with  abnormal  DTI  FA  measures  is  correlated  with  patients’ 

19  neurocognitive  status,  including  memory;  and  d)  DTI  FA  values  could  both  increase  and 

20  decrease  at  the  acute  setting,  which  is  also  a  new  finding  in  mTBI  research.  In  the  acute 

21  setting,  the  immediate  challenge  for  emergency  physicians  is  to  identify  those  CT 

22  negative  but  symptomatic  patients  with  intracranial  abnormalities  that  may  be  predictive 

23  of  long  term  neurocognitive  sequalae  [12].  Given  the  fact  that  most  mTBI  patients  stay  in 

24  the  emergency  department  for  a  few  hours  up  to  24  hours,  our  comprehensive  approach 

25  at  this  stage,  while  the  patients  are  still  in  the  emergency  department,  are  more  likely  to 
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1  help  emergency  physicians  make  decisions  on  patient’s  management  than  examining 

2  them  days  or  weeks  later. 

3 

4  This  study  ideally  extends  the  previous  work  by  our  co-authors  demonstrating  the 

5  relationships  between  different  pathways  for  UCH-L1  and  GFAP  and  different  types  of 

6  brain  injury  pathophysiology  after  severe  TBI  as  characterized  by  CT  [53].  In  addition  to 

7  the  previous  findings,  these  pilot  data  suggest  that  the  combined  use  of  biochemical 

8  markers  and  advanced  MRI  techniques  may  provide  an  important  tool  to  evaluate  and 

9  characterize  mTBI  patients  of  importance  for  the  understanding  of  the  different 

10  pathophysiological  mechanisms  following  TBI  and  for  the  development  of  effective 

11  therapies. 

12 

13  The  heterogeneity  of  brain  injury  pathology.  It  is  well  recognized  that  brain  injury 

14  pathology  is  heterogeneous  and  complex  [54],  Each  technique  employed  in  this  study 

15  brings  unique  aspects  of  brain  injury  pathology  and  contributes  to  the  whole  picture: 

16  intracranial  bleeding,  detected  by  SWI,  manifests  blood  vessel  damage  [25];  DTI  finding 

17  signifies  the  damage  of  WM  integrity  [13,  14];  UCH-L1  for  neuronal  injury  [44,  55];  and 

18  GFAP  for  glial  damage  [53,  56].  These  different  pathologies  may  be  correlated  with  each 

19  other  and,  together,  they  can  cover  the  spectrum  of  brain  injury  that  contributes  to 

20  impaired  brain  function.  Our  data  demonstrated  that  intracranial  bleeding  was  associated 

21  with  elevated  GFAP  levels,  in  suggestion  of  glial  injury  in  association  with  vascular 

22  damage.  Meanwhile,  the  non-association  between  UCPI-L1  and  MRI  data  implies  that 

23  neuronal  injury  may  not  happen  together  with  vascular  damage  or  axonal  injury. 

24  However,  they  all  demonstrated  abnormalities  of  mTBI  in  different  aspects,  in  suggestion 

25  that  they  are  also  complementary  to  each  other  for  brain  injury  detection.  This  further 

26  confirms  the  heterogeneity  of  brain  injury  pathology. 
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2  Intracranial  bleeding  and  elevated  GFAP  levels.  Searching  of  intracranial  bleeding  is 

3  critically  important  in  diagnostic  radiology.  The  confirmation  of  bleeding  in  mild  TBI  will 

4  automatically  categorize  a  patient  into  “complicated  mild  TBI”,  who  tends  to  have  worse 

5  outcome  than  those  without  any  intracranial  bleeding  [57],  A  most  recent  study  of  135 

6  mTBI  patients,  scanned  at  12  days  after  injury,  demonstrated  that  one  or  more  brain 

7  contusions  on  MRI,  and  >4  foci  of  hemorrhagic  axonal  injury  on  MRI,  were  each 

8  independently  associated  with  poorer  3-month  outcome  [58].  From  pathophysiological 

9  perspective,  GFAP  is  a  structural  protein  of  astroglial  cells  that  are  located  in  the 

10  intracellular  space  of  astrocytes.  The  damage  to  astrocytes  will  cause  the  release  of 

11  GFAP  into  extra-cellular  space  and  that  might  leak  into  the  blood  stream  through  a 

12  compromised  blood-brain  barrier  (BBB)  [59].  Furthermore,  the  end  processes  of 

13  astrocytes  surround  the  endothelial  cells  of  vasculature  system  and  make  astrocytes  an 

14  integral  part  of  neural  vascular  unit  [60].  The  damage  or  temporal  opening  of  the  BBB 

15  will  also  likely  further  damage  to  the  surrounding  astrocytes  as  well.  Supporting  this,  in 

16  stroke  studies,  considerable  amount  of  data  demonstrated  significantly  increase  of 

17  GFAP  in  expanding  intra-cerebral  hemorrhage  (ICH)  than  that  in  ischemic  stroke  [61, 

18  62],  Other  studies  even  reported  a  close  correlation  between  GFAP  serum  concentration 

19  and  ICH  volume  [63,  64],  Even  a  multi-center  clinical  trial  was  conducted  to  use  GFAP  to 

20  differential  ICH  from  ischemic  stroke  [63].  Instead  of  just  ICH,  our  data  demonstrated 

21  that  all  intracranial  hemorrhage  cases,  including  both  extra-axial  and  parenchymal 

22  hemorrhage,  have  significantly  elevated  GFAP  levels.  This  implies  that  GFAP  level  in 

23  blood  serum  has  the  potential  to  serve  as  a  quick  screening  biomarker  to  triage  mTBI 

24  patients  for  MRI  confirmation  of  intracranial  bleeding  for  an  unfavorable  outcome 

25  prediction. 

26 
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1  The  role  of  BBB.  The  elevated  biomarker  levels  measured  in  our  patients  support  the 

2  idea  of  a  BBB  breakdown  that  has  often  been  documented  in  patients  with  TBI  even 

3  after  mild  injuries  [65].  Indeed,  both  UCH-L1  and  GFAP  are  CNS-specific  proteins  with 

4  very  low  concentrations  in  blood  in  healthy  people,  almost  below  the  threshold  of 

5  detection  by  using  current  biomarker  technology  [40,  46].  The  elevated  level  of  either 

6  one  requires  the  same  pathway  to  leak  into  the  blood  pool:  the  breakdown  of  or 

7  compromised  BBB.  Given  the  much  smaller  size  of  UHCL-1  and  GFAP  than  red  blood 

8  cells,  these  proteins  could  more  easily  get  into  the  blood  stream  through  BBB  temporal 

9  opening.  Consequently,  the  elevated  biomarker  levels  seem  to  detect  a  BBB 

10  compromise  more  relevant  than  the  MRI-detectable  bleeding.  At  micro  level,  the  BBB 

11  damage  may  not  be  severe  enough  or  the  temporal  opening  of  BBB  may  not  be  long 

12  enough  to  allow  enough  red  blood  cells  to  pass  through  or  cause  sufficient  amount  of 

13  leakage  that  makes  it  visible  as  bleeding  on  neuroimaging  at  the  macro-level  [40]. 

14  Flowever,  this  BBB  compromise  may  be  already  big  enough  for  sufficient  amount  of 

15  small  protein  biomarkers  leaking  into  the  blood  stream  and  make  it  detectable  in  modern 

16  biomarker  technique.  Compared  with  detectable  bleeding,  which  consists  of  only  a  small 

17  fraction  of  mTBI  patients,  the  elevation  of  CNS-specific  proteins  in  the  blood  pool  might 

18  be  able  to  serve  as  a  more  sensitive  biomarker  for  the  compromise  of  BBB  in  mTBI  at 

19  the  acute  stage. 

20 

21  Correlation  with  patients’  neurocognitive  performance.  Our  data  showed  that,  DTI 

22  lesion  loads,  both  TBSS  FA  and  VBA  FA  lesion  loads,  are  correlated  with  their  SAC 

23  score  and  delayed  recall.  More  evidence  reported  that  DTI  FA  values  are  correlated  with 

24  mTBI  patients’  neurocognitive  outcome  [14,  66].  Particularly,  certain  regions  of  DTI  WM 

25  tract  are  correlated  with  patients’  specific  neurocognitive  outcome  [14,  66].  As  an 

26  example,  Muhkerjee  et  al  reported  that  DTI  findings  are  correlated  with  patients’ 
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1  neurocognitive  performance,  but  not  hemorrhage  [67,  68].  Our  DTI  finding  at  the  acute 

2  stage  is  in  the  same  line  as  the  published  result  at  the  sub-acute  or  chronic  stage.  This 

3  further  confirms  the  hypothesis  that  there  might  be  microstructural  damages  or  changes 

4  in  WM  tracts  that  account  for  patients’  neurocognitive  deficits.  However,  this  small-scale 

5  damage  may  not  reach  to  the  degree  of  vessel  rupture  yet  that  causes  bleeding  or 

6  hemorrhage. 

7 

8  DTI  FA  increase  or  decrease.  Two  more  DTI  findings  in  this  study  are  new  to  the  field: 

9  a)  the  co-existence  of  both  increased  and  decreased  FA  values  in  mTBI  patients  within 

10  24  hours  after  injury,  and  b)  the  dominance  of  increased  FA  lesions.  All  DTI  studies  of 

11  moderate  to  severe  TBI  patients  [23,  69-71]  and  subacute/chronic  mTBI  patients  [22,  67, 

12  72-74]  report  FA  decreases  which  are  correlated  with  clinical  or  neuropsychological 

13  measures.  However,  there  are  seemingly  contradictory  findings  in  mild  TBI  in  the  acute 

14  stage  (within  one  week  after  injury)  in  the  literature:  Inglese  [22]  and  Arfanakis[21]  both 

15  reported  FA  decreases,  while  Wilde  [75],  Bazarian  [76],  and  Mayer  [77]  reported  FA 

16  increases  and  decreased  radial  diffusivity.  Furthermore,  Michael  Lipton  et  al  [78] 

17  reported  bi-directional  changes  (both  increase  and  decrease)  of  FA  in  chronic  mTBI 

18  patients.  Of  particular  note,  the  terminology  “acute  stage”  could  mean  quite  different 

19  timing  frames  across  the  studies:  some  defined  it  as  within  24  hours  after  injury  and 

20  some  even  as  within  7  days  after  injury.  To  date,  only  two  studies  reported  MRI  scan  of 

21  mTBI  patients  within  24  hours  after  injury  and  they  all  have  only  a  handful  of  patients 

22  [21 , 75].  The  co-existence  of  both  FA  decrease  and  increase  within  24  hours  after  injury 

23  is  a  new  finding  in  the  field.  It  further  demonstrates  the  heterogeneity  of  mTBI  pathology 

24  at  this  stage.  Meanwhile,  lesion  load  with  increased  FA  is  significantly  higher  than  that 

25  with  decreased  FA.  Increased  FA  lesion  load,  but  not  decreased  FA  lesion  load,  is 

26  correlated  with  patients’  SAC  and  delayed  recall  scores  in  our  data.  It  has  been 
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1  suggested  that  increased  FA  acutely  may  reflect  cytotoxic  oedema  [75],  which  would 

2  shunt  extracellular  fluid  into  swollen  cells.  This  could  have  the  effect  of  reducing  inter- 

3  axonal  free  water  and  therefore  increasing  anisotropy.  This  demonstrated  that  cytotoxic 

4  edema  might  be  the  major  pathology  that  accounts  for  patients’  neurocognitive 

5  symptoms  at  hyper-acute  stage  (within  24  hours  after  injury). 

6 

7  The  need  for  an  axonal  injury  biomarker.  Our  findings  confirm  that  a  panel  of 

8  biomarkers  rather  than  a  single  analyte  seem  to  have  the  most  utility  for  the  diagnosis  of 

9  mTBI  patients,  and  improved  characterization  of  the  injury.  Importantly,  in  the  current 

10  study  neither  UCH-L1  nor  GFAP  was  associated  with  WM  injury  identified  by  DTI.  Since 

11  the  traumatic  axonal  injury  is  believed  to  be  a  major  determinant  of  functional  and 

12  neurocognitive  symptoms  following  TBI  as  demonstrated  by  the  correlation  between  DTI 

13  and  patients’  neurocognitive  deficits,  there  might  be  a  need  for  specific  axonal  injury 

14  biomarkers.  Further  work  is  needed  to  develop  additional  biomarker  platforms,  including 

15  axonal  injury  markers,  in  addition  to  the  neuronal  and  glial  damage  proteins  examined 

16  here,  and  to  identify  the  relationships  with  advanced  MRI  techniques  and  patient 

17  outcomes  that  will  help  validate  and  confirm  their  clinical  utilities  in  the  acute  setting  [1 2], 

18 

19  Limitations  and  future  work.  Despite  its  encouraging  finding,  this  preliminary  work  has 

20  limitations,  including  a  small  sample  size  and  the  lack  of  long-term  outcome  data. 

21  Additional  research  will  be  required  to  validate  our  current  findings  in  a  large  cohort  of 

22  patients  with  longitudinal  follow  up  and  to  further  determine  the  relationships  between 

23  neuroimaging  and  biomarker  findings  in  the  prediction  of  mTBI  outcome. 

24 

25  CONCLUSIONS 
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1  To  summarize,  this  work  represents  the  first  effort  of  combining  both  blood  protein 

2  biomarkers  and  advanced  MRI  to  improve  the  detection  and  characterization  of  brain 

3  injuries  after  mild  TBI  in  the  acute  stage  (within  24  hours  after  injury).  Our  data 

4  demonstrate  elevated  GFAP  and  UCH-L1  levels  in  mTBI  patients  at  the  acute  stage  in 

5  comparison  with  controls.  Particularly,  all  cases  with  intra-cranial  hemorrhage  had 

6  significantly  higher  GFAP  levels  than  those  without  hemorrhage.  Patients’  DTI  measures 

7  were  correlated  with  their  neurocognitive  status  at  this  stage.  This  overlapping  and 

8  complementary  role  of  blood  biomarkers  and  imaging  to  brain  injury  detection  offers  the 

9  promise  that  they  might  be  used  in  conjunction  in  the  management  of  patients  with 

10  mTBI.  Further  studies  with  larger  numbers  of  patients  will  be  required  to  assess  the 

11  reproducibility  of  these  findings  and  to  confirm  the  potential  clinical  utilities  as  diagnostic 

12  adjuncts  in  the  acute  setting. 

13 

14 
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1 


FIGURE  LEGENDS 


2  Figure  1.  Correlations  between  DTI  lesion  load  and  patients’  neurocognitive  data. 

3  As  demonstrated  in  the  figures,  DTI  lesion  load  (both  TBSS  and  VBA  data)  are 

4  significantly  correlated  with  patients’  overall  SAC  score  and  delayed  recall.  R  squared 

5  values  are  shown  on  each  figure  for  linear  regression. 

6 

7  Figure  2.  Dots  plots  demonstrating  UCH-L1  and  GFAP  concentrations.  Serum 

8  UCH-L1  (A)  and  GFAP  (B)  concentrations  on  admission  in  TBI  patients  and  in  controls. 

9  Error  bars  represent  median  and  range.  Significant  differences  are  indicated  with  **  (P< 

10  0.01 )  or  ***  (P  <0.001 )  (Mann-Whitney  U-test) 

11 

12  Figure  3.  Box-and-whisker  plots  demonstrating  UCH-L1  and  GFAP  concentrations. 

13  (A)  Serum  UCH-L1  and  GFAP  concentrations  in  patients  who  were  victims  of  assault 

14  and  in  patients  injured  in  a  MVA.  (B)  Serum  UCFI-L1  and  GFAP  concentrations  in 

15  patients  with  ventricular  hemorrhages  and  hemorrhagic  contusions  and  in  patients  with 

16  non-hemorrhagic  lesions.  The  black  horizontal  line  in  each  box  represents  the  median, 

17  with  the  boxes  representing  the  interquartile  range.  Significant  differences  are  indicated  * 

18  (P<  0.05)  or  **  (P  <0.01)  (Mann-Whitney  U-test). 
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1  Figure  4.  Case  1 .  MRI  and  Biomarker  Profile  in  a  Patient  with  Intraventricular 

2  Hemorrhage  but  Missed  by  CT.  Panels  a)  and  b)  are  SWI  images  at  different  locations 

3  of  the  brain  showing  intra-ventricular  blood  and  left  lingual  gyrus  blood  product  (see 

4  arrows);  panel  c)  is  FLAIR  image  showing  non-specific  white  matter  hyper-  intensities 

5  (see  arrows);  panel  d)  is  DTI  FA  map  showing  the  co-existence  of  voxels  with  increased 

6  and  decreased  FA  measures  (red  color  means  FA  decrease  and  blue  color  FA  increase 

7  in  comparison  with  controls,  t>3  for  t-test);  and  panel  e)  is  blood  biomarker  temporal 

8  profile,  which  exhibiting  extraordinarily  high  GFAP  levels  over  time  in  comparison  with 

9  controls  (median  0.004,  interquartile  range  0.004-0.015).  Despite  being  missed  by  CT, 

10  the  patient  case  was  still  detected  by  both  blood  biomarker  and  MRI. 

11 

12  Figure  5.  Case  2.  MRI  and  Biomarker  Profile  in  a  Patient  with  traumatic  axonal 

13  injury  but  normal  appearing  structural  MRI.  Panels  a-d)  are  MRI  images  at  the 

14  corpus  callosum  and  fornix  level.  Panels  e-h)  are  MRI  images  at  the  level  of  superior 

15  coronal  radiata.  Panel  i)  is  blood  biomarker  profile.  FLAIR  and  SWI  images  both  indicate 

16  the  skulp  contusion  at  the  parieto-occipital  region  (long  arrows)  but  normal  appearing 

17  brain  structure.  Flowever,  both  DTI  TBSS  and  VBA  analyses  detected  significantly 

18  reduced  FA  values  at  the  ipsilateral  side  (corticospinal  tract)  and  contralateral  side 

19  (superior  corona  radiata)  of  brain  white  matter  (arrow  heads),  in  suggestion  of  coup  and 

20  contra  coup  injury  at  the  microstructure  of  white  matter.  Cold  color  indicates  reduced  FA 

21  values  in  comparison  with  controls  (t>3  for  t-test).  Blood  biomarkers  indicate  slightly 

22  increased  GFAP  levels  over  time  but  significantly  increased  UCPI-L1  at  the  admission. 

23 

24  Figure  SI.  Serum  biomarker  levels  over  the  first  24  hours  after  mTBI  compared 

25  with  controls.  Serum  UCPI-L1  (A)  levels  are  maximal  early  after  injury  (on  admission) 
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1  (TE=0.24  [0.096-0.346]),  while  GFAP  (B)  concentrations  peaked  12  hours  after  injury 

2  (0.35  [0.036-2.56]).  Error  bars  represent  median  and  IQR. 

3 
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1 

2  Table  1 .  Individual  demographic  and  clinical  data  of  the  9  patients  enrolled  in  the 

3  study 


Patient 

no. 

Age/ 

Gender 

Race 

GCS 

Mechanis 
m  of  Injury 

SA 

C 

CT 

Structural  MRI 

P-001 

56/M 

Asian 

15 

MVA 

Negative 

Nonspecific  WM 
hyperintensies,  Small  foci 
of  intraventricular  blood 
on  the  left.  Small  blood 
product  in  the  left  lingual 
gyrus 

P-002 

36/F 

Black 

15 

MVA 

23 

Negative 

Negative 

P-003 

19/M 

Black 

15 

MVA 

25 

Negative 

Nonspecific  FLAIR 
hyerintensity  in  posterior 
cerebral  WM 

P-004 

35/M 

Arabic 

15 

Assault 

26 

Positive 
hemorrhagic 
contusion  on 
left  parieto¬ 
temporal  lobe 

Hemorrhagic  contusion 
on  left  parieto-temporal 
lobe,  left  ventricular 
hemorrhage 

P-005 

52/M 

Black 

15 

Assault 

19 

Negative 

Non-specific,  multiple 
scattered  discrete  foci  in 
cerebral  WM 

P-006 

53/M 

Black 

15 

MVA 

24 

Negative 

Non-specific:  Two 
isolated  punctate  foci  of 
blood  in  the  right  periatrial 
WM 

P-007 

39/M 

Cauca 

sian 

15 

assault 

22 

Small  SAH  in 
right  Sylvian 
Fissure 

SAH  in  right  Sylvian 

Fissure 

P-008 

58/M 

Cauca 

sian 

15 

Assault 

19 

Negative 

Non-specific,  super  sella 
lesion,  congenital  cistern 
lesion  in  posterial  fossa 

P-009 

23/M 

Cauca 

sian 

15 

Assault 

24 

Negative 

Negative 

10 
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1 

2  Table  2.  Serum  concentration  of  UCH-L1  and  GFAP  in  patients  with  mTBI  and  in 

3  controls.  Data  are  given  as  median  (interquartile  range). 

4 


Serum  UCH-L1  (ng/mL) 

Serum  GFAP  (ng/mL) 

TBI 

Admission 

0.242  (0.096-0.336)* 

0.043  (0.015-0.375)* 

Hemorrhagic 

0.164(0.098-0.314)* 

0.517(0.239-4.610)**  f 

Non-Hemorrhagic 

0.171(0.107-0.248)** 

0.015(0.015-0.06)** 

Controls 

0.05  (0.05-0.153) 

0.004  (0.004-0.015) 

5  *  p  <.01  and  **  p<.001  (p  values  of  the  Mann-Whitney  test  for  differences  between  the  groups  [TBI  versus 

6  Controls]) 

7  t  p<-01  (p  values  of  the  Mann-Whitney  test  for  differences  between  the  groups  [Hemorrhagic  versus  Non- 

8  Hemorrhagic]) 

9 

10 

11 

12 

13 

14 

15 
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Figure  1 

Click  here  to  download  high  resolution  image 
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Quantitative  Susceptibility  Mapping  of  Small  Objects 
Using  Volume  Constraints 

Saifeng  Liu,1  Jaladhar  Neelavalli,2*  Yu-Chung  N.  Cheng,2  Jin  Tang,1 
and  E.  Mark  Haacke1-3 


Microbleeds  have  been  implicated  to  play  a  role  in  many  neu¬ 
rovascular  and  neurodegenerative  diseases.  The  diameter  of 
each  microbleed  has  been  used  previously  as  a  possible  quan¬ 
titative  measure  for  grading  microbleeds.  We  propose  that 
magnetic  susceptibility  provides  a  new  quantitative  measure  of 
extravasated  blood.  Recently,  a  Fourier-based  method  has 
been  used  that  allows  susceptibility  quantification  from  phase 
images  for  any  arbitrarily  shaped  structures.  However,  when 
very  small  objects,  such  as  microbleeds,  are  considered,  the 
accuracy  of  this  susceptibility  mapping  method  still  remains  to 
be  evaluated.  In  this  article,  air  bubbles  and  glass  beads  are 
taken  as  microbleed  surrogates  to  evaluate  the  quantitative  ac¬ 
curacy  of  the  susceptibility  mapping  method.  We  show  that 
when  an  object  occupies  only  a  few  voxels,  an  estimate  of  the 
true  volume  of  the  object  is  necessary  for  accurate  susceptibil¬ 
ity  quantification.  Remnant  errors  in  the  quantified  susceptibil¬ 
ities  and  their  sources  are  evaluated.  We  show  that  quantifying 
magnetic  moment,  rather  than  the  susceptibility  of  these  small 
structures,  may  be  a  better  and  more  robust  alternative.  Magn 
Reson  Med  69:716-723,  2013.  ©2012  Wiley  Periodicals,  Inc. 

Key  words:  susceptibility  mapping;  microbleeds;  air  bubbles; 
magnetic  moment 


The  measurement  of  magnetic  susceptibility  offers  an 
entirely  new  form  of  contrast  in  magnetic  resonance 
imaging  (1-6).  More  specifically,  susceptibility  quantifi¬ 
cation  has  already  found  applications  in  mapping  out 
iron  in  the  form  of  ferritin  in  brain  tissues  such  as  the 
basal  ganglia  (1,4,6)  and  in  the  form  of  deoxyhemoglobin 
for  measuring  the  oxygen  saturation  in  veins  (1,4).  This 
new  form  of  imaging  may  provide  a  means  for  monitor¬ 
ing  longitudinal  changes  in  iron  content  in  dementia, 
multiple  sclerosis,  traumatic  brain  injury,  and  Parkin¬ 
son’s  disease.  It  may  also  be  used  to  monitor  microbleeds 
which  have  been  implicated  in  the  progression  of  vascu¬ 
lar  dementia  (7),  Alzheimer’s,  and  other  neurovascular 
disorders  (8,9). 
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One  of  the  most  recent  susceptibility  mapping  meth¬ 
ods  is  a  Fourier-based  method  (2,3,5,10)  which  utilizes 
phase  images.  The  accuracy  of  such  a  method  depends 
on  the  volume  measurement  of  the  object.  For  example, 
to  quantify  the  susceptibility  of  a  given  microhleed,  usu¬ 
ally  the  center  and  the  radius  of  the  microbleed  have  to 
be  determined  (10-13).  Alternate  volume  estimations  of 
the  microbleed  from  high-resolution  spin  echo  images 
may  overcome  these  limitations.  With  a  gradient  echo 
sequence,  the  apparent  volume  of  the  object  is  increased 
owing  to  what  is  commonly  referred  to  as  the  “bloom¬ 
ing”  effect,  a  signal  loss  around  the  object  caused  by  T2* 
dephasing.  This  increased  apparent  volume  may  be  used 
to  obtain  an  estimate  of  susceptibility  while  the  product 
of  the  apparent  volume  and  the  estimated  susceptibility 
is  much  more  robust  and  should  still  provide  a  good 
estimate  of  the  magnetic  moment  of  the  object. 

The  goal  of  this  article  is  to  evaluate  the  quantitative 
accuracy  of  a  Fourier-based  susceptibility  mapping 
method  when  it  is  applied  to  small  structures,  and  to 
show  that:  (1)  an  accurate  estimate  of  the  magnetic 
moment  is  possible  using  multiecho  gradient  echo  imag¬ 
ing;  and  (2)  the  accuracy  of  the  effective  susceptibility 
can  be  improved  using  the  magnetic  moment  when  an 
estimate  of  the  true  volume  is  available.  For  validation, 
we  used  a  gel  phantom  with  air  bubbles  and  glass  beads 
to  mimic  the  clinical  situation  of  microbleeds.  The 
method  illustrated  here  does  not  depend  on  the  suscepti¬ 
bility  value  or  the  size  of  the  object. 


THEORY 

Current  susceptibility  mapping  methods  are  based  on  the 
relationship  between  the  susceptibility  distribution  and  the 
magnetic  field  variation  in  the  Fourier  domain  (1-6,10): 

A  B(k)  =  G(k)Ax(k)  [1] 

where  A B(k)  is  the  Fourier  transform  of  the  magnetic 
field  variation  A B(r),  A x(k)  is  the  Fourier  transform  of  the 
susceptibility  distribution  Ax(r),  and  G(k )  is  the  Green’s 
function 

G(k)  =  1/3 -k2/(k2  +  tf+k2J  [2] 

assuming  that  the  main  field  direction  is  in  the  z-direc- 
tion.  Susceptibility  quantification  is  an  ill-posed  inverse 
problem,  owing  to  zeros  in  the  Green’s  function  G(k) 
along  the  magic  angle  in  the  k-space  domain.  As  a  result, 
regularization  is  required,  as  G(k)_1  is  needed  in  the 
evaluation  of  susceptibility.  In  this  study,  we  applied  the 
regularization  procedure  described  in  a  previous  study 
(1)  in  which  the  intensity  of  G(k)-1  is  reasonably  attenu¬ 
ated  when  the  absolute  value  of  G(k)  is  below  a 
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threshold  value.  The  selection  of  this  threshold  is  a 
trade-off  between  the  susceptibility-to-noise  ratio  of  the 
reconstructed  susceptibility  map  and  the  accuracy  in 
susceptibility  quantification  (1,6).  The  threshold  value 
was  chosen  to  be  0.1  in  this  study. 

Although  ideally  Ax  is  the  sought  after  parameter, 
when  reduced  resolution  or  T2*  effects  confound  a  clean 
measurement  of  the  object’s  volume,  it  is  more  appropri¬ 
ate  to  investigate  the  associated  magnetic  moment  (or, 
equivalently,  the  total  or  integrated  susceptibility 
weighted  by  the  voxel  volume  (12)).  To  see  why  this  is 
the  case,  consider  a  sphere  with  a  susceptibility  differ¬ 
ence  Ax,  the  induced  magnetic  field  at  point  P(r,  0)  out¬ 
side  the  sphere  is  given  by  (14): 


aP  M  AXr3(3cos2e-l)B0 
AAout(r)  - - —5 - 


[3] 


where  Ax  =  Xin  —  Xout,  Xin  is  the  susceptibility  inside  the 
object,  Xout  is  the  susceptibility  outside  the  object,  r0  is 
the  radius  of  the  sphere,  r  is  the  distance  from  the  point 
P(r,  0)  to  the  center  of  the  sphere,  and  0  is  the  angle 
between  the  point  P  and  the  main  field  direction.  For 
simplicity,  the  meaning  of  susceptibility  in  this  article 
will  be  taken  to  be  Ax  rather  than  Xin  or  xout-  Equation  3 
also  indicates  that  the  product  AyR  is  independent  of 
echo  time  (TE),  where  V  =  4m^/3  is  the  true  volume  of 
the  sphere.  The  magnetic  dipole  moment  of  the  spherical 
object  is  given  as  (14): 

=  4tt t30M  _  4TTrgAX50 
^  3  ~  (3p0)  [  1 

when  Ax  is  much  smaller  than  1.  Here  jjl0  is  the  perme¬ 
ability  of  free  space  and  M  is  the  induced  magnetization. 
The  product  of  Ax  and  volume  V  is  the  effective  magnetic 
moment  term  (AxV)  and  is  simply  referred  to  as  magnetic 
moment  hereafter  in  this  paper.  The  phase  value  at  a  par¬ 
ticular  TE  is  given  in  a  right-handed  system  by: 

A<j>(r)  = -yABout(r)TE  [5] 

The  susceptibility  Ax  may  be  quantified  using  the  phase 
information  if  the  true  volume  (V)  of  the  object  is  known. 
Otherwise,  the  magnetic  moment  (AxV)  may  be  found.  As 
gradient  echo  images  lead  to  a  dephasing  artifact  and  the 
object  appears  larger  than  its  actual  size,  we  defined  an 
apparent  volume  V  and  assuming  that  the  susceptibility 
of  this  larger  object  can  be  accurately  quantified,  the  mag¬ 
netic  moment  could  still  be  accurately  calculated.  An  esti¬ 
mated  susceptibility  value  Ax'  can  be  calculated  from  the 
Fourier-based  method  using  Eqs.  1-5.  The  quantity  Ax'V' 
provides  an  estimate  of  the  magnetic  moment.  Finally,  the 
true  susceptibility  Ay  can  be  calculated  using  the  follow¬ 
ing  equation: 


AX  =  Ax'V'/V  [6] 

In  this  study,  we  use  three  volume  definitions.  The 
first  one  is  the  true  volume  V.  The  second  one  is 
the  apparent  volume  V ,  which  is  used  in  estimating  the 


magnetic  moment.  The  apparent  volume  is  related  to  the 
signal  loss  owing  to  T2*  dephasing  and  is  determined 
from  gradient  echo  magnitude  images,  as  described  later. 
The  last  one  is  the  spin  echo  volume  Vse,  which  is  meas¬ 
ured  from  the  spin  echo  images.  This  volume  is  used  as 
an  MR-based  estimate  of  the  true  volume.  We  used  simu¬ 
lations  and  multiecho  gradient  echo  images  of  a  gel 
phantom  containing  air  bubbles  and  glass  beads  of 
varying  sizes  to  test  Eq.  6.  Although  glass  beads  can  be 
considered  as  almost  perfect  spheres,  air  bubbles  are 
closer  to  the  clinical  situation  of  variable-shaped 
microbleeds. 

METHODS 

Simulations 

To  evaluate  the  validity  of  Eq.  6  for  susceptibility  calcu¬ 
lation  of  small  objects,  we  simulated  magnitude  and 
phase  images  of  four  spheres  with  different  radii  at  21 
different  TEs  (from  0  to  20  ms,  with  a  step  size  of  1  ms). 
In  each  simulation,  the  sphere  was  placed  in  the  center 
of  a  10243  matrix  with  complex  elements.  The  radii  of 
four  spheres  tested,  within  this  10243  matrix,  were  32, 
48,  64,  and  96  pixels,  respectively.  The  magnitude  inside 
each  sphere  was  set  to  0,  whereas  the  background  magni¬ 
tude  was  set  to  300  to  simulate  intensities  in  the  experi¬ 
mental  data  from  the  gel  phantom.  The  phase  images  of 
the  spheres  were  generated  according  to  Eqs.  3  and  5 
with  Ax  =  9.4  ppm.  To  simulate  Gibbs  ringing  as  well  as 
partial  volume  effects  seen  in  actual  MR  data,  a  process 
simulating  the  MR  data  sampling  was  used.  Complex 
images  generated  in  each  10243  matrix  were  Fourier 
transformed  into  k-space.  The  central  323  region  was 
selected  from  k-space  and  was  inverse  Fourier  trans¬ 
formed  back  to  the  imaging  domain  generating  low-reso¬ 
lution  data  containing  both  Gibbs  ringing  and  partial-vol¬ 
ume  effects.  The  radii  of  the  four  spheres  became  1,  1.5, 
2,  and  3  pixels,  respectively,  in  this  final  32s  volume. 
White  gaussian  noise  was  then  added  to  the  real  and 
imaginary  channels  of  the  complex  data  in  the  image  do¬ 
main  such  that  the  signal-to-noise  ratio  in  resultant  mag¬ 
nitude  images  was  10:1.  Susceptibility  and  the  magnetic 
moment  values  were  quantified  for  each  of  the  spheres 
at  all  TEs  and  errors  associated  with  these  measurements 
were  evaluated. 

Phantom  Experiments 

A  gel  phantom,  containing  14  small  air  bubbles  and  9 
glass  beads  of  varying  sizes,  was  imaged  at  3T  (Siemens 
VERIO,  Erlangen,  Germany)  using  a  five-echo  3D  gradi¬ 
ent  echo  sequence.  The  TEs  were  3.93,  9.60,  15.27, 
20.94,  and  26.61  ms,  respectively.  Other  imaging  parame¬ 
ters  for  the  gradient  echo  sequence  were  repetition  time 
33  ms,  flip  angle  11°,  read  bandwidth  465  Hz/pixel, 
voxel  size  0.5  x  0.5  x  0.5mm3,  and  matrix  size  512  x 
304  x  176.  A  multislice  2D  spin  echo  data  set  was  also 
collected  with  flip  angle  of  90°,  repetition  time  of  5000 
ms,  and  TE  =  15  ms  and  with  the  same  field  of  view, 
bandwidth,  resolution,  and  matrix  size  as  in  the  gradient 
echo  data  set.  This  is  to  maintain  a  one-to-one  correspon¬ 
dence  of  the  spin  echo  with  the  gradient  echo  images  of 
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the  phantom.  To  ensure  that  the  field  perturbation  meas¬ 
ured  in  the  phase  images  is  the  actual  perturbation  pro¬ 
file  from  the  gel  phantom,  we  first  performed  shimming 
using  a  spherical  phantom  immediately  before  perform¬ 
ing  the  imaging  experiment.  Manual  shimming  was  per¬ 
formed  on  the  spherical  phantom,  to  a  spectral  full- 
width  at  half-maximum  of  13  Hz  and  the  shim  coeffi¬ 
cients  were  noted.  The  same  shim  settings  were  used 
while  imaging  the  gel  phantom  to  ensure  that  field  per¬ 
turbation  profile  owing  to  the  presence  of  the  phantom 
in  the  magnet  is  not  influenced  by  any  additional 
shimming. 

For  the  construction  of  the  phantom,  an  agarose  gel  so¬ 
lution  was  prepared  with  an  8%  concentration  by  weight 
and  poured  into  a  cylindrical  container.  In  the  lower 
portion  of  the  container,  the  gel  was  first  filled  to  one- 
third  the  height  of  the  cylinder  and  nine  glass  beads  of 
various  sizes  were  embedded  in  the  gel.  The  true  diame¬ 
ter  of  the  glass  beads  was  roughly  measured  using  cali¬ 
pers  before  the  glass  beads  were  put  into  the  gel  solu¬ 
tion.  Specifically,  four  glass  beads  were  2  mm  in 
diameter,  three  glass  beads  were  3  mm  in  diameter,  one 
glass  bead  was  5  mm  in  diameter,  and  the  largest  glass 
bead  was  6  mm  in  diameter.  The  phantom  was  allowed 
to  cool  so  that  the  gel  solidified  and  properly  engulfed 
the  glass  beads.  The  rest  of  the  prepared  gel  solution  was 
then  poured  into  the  cylindrical  container  and  variable¬ 
sized  bubbles  were  injected  by  pumping  various  amount 
of  air  into  the  gel  using  an  empty  syringe  (two  smallest 
air  bubbles  were  excluded  from  this  study,  owing  to  the 
limitation  in  volume  estimation  of  small  objects;  details 
are  provided  in  later  sections).  The  theoretical  suscepti¬ 
bility  difference  between  air  and  water  is  known  to  be 
9.4  ppm  and  will  be  used  to  compare  with  the  measure¬ 
ments  from  our  method.  For  glass  beads,  the  susceptibil¬ 
ity  values  were  measured  independently  in  a  former 
study  to  be  —1.8  ±  0.3  ppm  relative  to  water  (15). 

First,  to  identify  air  bubbles  and  glass  beads  in  the  col¬ 
lected  MR  data,  binary  masks  from  magnitude  data  were 
used.  The  intensity  variation  in  the  magnitude  images 
caused  by  the  RF  field  inhomogeneity  was  first  removed 
using  a  2D  quadratic  fitting,  before  the  binary  masks 
were  created.  A  reasonably  uniform  magnitude  intensity 
profile  across  the  phantom  was  obtained  after  this  inten¬ 
sity  correction.  The  binary  masks  were  created  by  local 
thresholding  of  the  corrected  magnitude  images  (11). 
First,  a  relatively  strict  threshold  was  used  to  pick  only 
the  voxels  where  the  signal  was  <50%  of  the  signal  in 
the  gel  away  from  the  air  bubbles  or  glass  beads,  as  both 
air  bubbles  and  glass  beads  have  much  lower  intensities 
than  the  intensity  of  the  surrounding  gel.  Next,  the  mean 
(amag-gei)  and  standard  deviation  (amag_gei)  were  calcu¬ 
lated  for  a  cubic  21  x  21  x  21  voxels  volume  of  interest 
for  each  bubble  or  glass  bead.  A  voxel  roughly  at  the  cen¬ 
ter  of  the  bubble  or  glass  bead  was  first  chosen  to  center 
this  213  voxel  window.  The  voxels  picked  up  in  the  first 
step  were  excluded  in  the  mean  and  standard  deviation 
calculation.  If  a  neighboring  voxel  has  intensity  lower 
than  amag_gei  —  Pcrmag_gei,  it  was  regarded  as  a  voxel 
belonging  to  air  or  glass  bead.  For  the  high  signal-to- 
noise  ratio  data  used  here,  (3  was  empirically  chosen  to 
be  4  to  separate  air  bubbles  and  glass  beads  from  gel. 


Susceptibility  Quantification 

To  reduce  the  background  field  or  phase  variation,  a  for¬ 
ward  modeling  approach  was  used  to  estimate  air/gel- 
phantom  interface  effects  (16).  The  phase  processing 
steps  were  as  follows: 

i.  The  original  phase  images  were  first  unwrapped 
using  the  phase  unwrapping  tool,  PRELUDE,  in 
FMRIB  Software  Library  (FSL)  (17).  With  the  geom¬ 
etry  of  the  gel  phantom  extracted  from  the  magni¬ 
tude  images  at  the  shortest  TE  (3.93  ms  in  this 
study),  the  background  field  effects  were  reduced 
by  fitting  the  predicted  phase  to  the  unwrapped 
phase  by  a  least  squares  method.  An  additional  2D 
quadratic  fitting  was  added  to  remove  the  induced 
phase  owing  to  eddy  currents. 

ii.  The  phase  value  inside  a  particular  air  bubble/glass 
bead  (where  the  binary  mask  is  1)  was  set  to  the 
mean  phase  (essentially  0)  from  the  local  21  x  21 
x  21  window  excluding  the  voxels  which  belong 
to  the  air  bubble  or  glass  bead.  This  is  owing  to  the 
fact  that  the  phase  inside  a  sphere  is  theoretically 
zero  and  the  nonzero  phase  is  induced  by  the  rem¬ 
nant  background  field  variation  as  well  as  Gibbs 
ringing.  This  step  also  determines  the  apparent 
volume  [V)  from  magnitude  images. 

iii.  At  each  TE,  a  160  x  160  x  87  voxel  volume  was 
cropped  from  the  original  phase  images.  This  vol¬ 
ume  was  selected  because  it  covers  most  of  the  air 
bubbles  and  glass  beads,  whereas  voxels  near  the 
edge  of  the  gel  phantom  were  excluded.  The 
selected  volume  was  then  zero-filled  to  a  512  x 
512  x  256  matrix. 

iv.  Susceptibility  maps  were  generated  using  a  thresh¬ 
old-based  approach  described  previously  in  Ref. 
(1).  The  mean  (ax.air  or  ax_giass)  and  standard  devia¬ 
tion  (CTx_air  or  <Tx.giass)  of  the  susceptibility  values  of 
air  bubble  (or  glass  bead)  were  measured,  taking 
into  account  the  background  susceptibility  of  the 
gel.  Measurements  were  obtained  in  the  following 
manner:  the  background  mean  (ax_gei)  and  standard 
deviation  (crx-gei)  of  the  local  gel  susceptibility 
value  around  each  bubble  or  glass  bead  was  first 
calculated  from  the  213  voxel  region  centered 
around  each  of  the  bubble/bead.  Within  this  213 
volume,  the  voxels  belonging  to  the  bubble  or 
glass  bead,  as  determined  by  the  binary  mask, 
were  excluded  for  this  background  mean  and 
standard  deviation  calculation.  Once  these  meas¬ 
ures  were  obtained,  for  susceptibility  of  air  bub¬ 
bles,  only  voxels  with  susceptibility  values  higher 
than  ax.gei  +  3-ax.gei  were  used  for  calculation  pur¬ 
poses;  whereas  for  glass  beads,  only  voxels  with  a 
susceptibility  value  lower  than  ax_gei  —  3crx_gei 
were  used.  This  process  assumes  that  the  noise  in 
the  susceptibility  maps  follows  a  gaussian  distri¬ 
bution,  and  the  susceptibility  of  a  voxel  consisting 
of  air  or  glass  is  statistically  different  from  a  voxel 
consisting  of  gel.  The  change  in  sign  is  owing  to 
the  fact  that  the  air  bubbles  are  paramagnetic  rela¬ 
tive  to  the  gel,  whereas  glass  beads  are  diamag¬ 
netic.  To  account  for  the  baseline  shift  caused  by 
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FIG.  1.  Apparent  volume  normalized  to  the  volume  at  TE  =  20  ms  (first  column),  measured  susceptibility  (second  column),  and  normal¬ 
ized  magnetic  moments  (third  column)  measured  at  different  TEs  of  four  different  spheres.  The  dashed  lines  in  the  second  column  (b,  e, 
h,  and  k)  indicate  the  true  susceptibility  9.4  ppm.  For  each  sphere,  the  effective  magnetic  moments  were  normalized  to  the  true  effec¬ 
tive  magnetic  moment. 


remnant  field  variation,  the  susceptibility  of  the 
air  bubble  (or  glass  bead)  was  taken  as  ax_air  — 

aX-gel  ax_glasS  ax-gel)- 

Volume  Measurement 

The  apparent  volume  of  the  air  bubble  or  glass  bead  was 
determined  from  the  binary  masks  directly,  i.e.,  by 
counting  the  number  of  voxels  inside  the  air  bubble  or 
glass  bead.  On  the  other  hand,  the  spin  echo  volume  is 


measured  utilizing  the  “object  strength”  notion  proposed 
by  Tofts  et  al.  (18),  in  which  the  total  intensity  is  meas¬ 
ured  for  a  particular  volume  of  interest.  For  a  volume 
composed  of  two  types  of  tissues,  a  and  b,  the  total  in¬ 
tensity  can  be  expressed  as 

I  =  7a  •  na  +  Ih  ■  (N  -  na)  =  (7a  -  Ib)  ■  na  +  Ib  ■  N  [7] 

where  I  is  the  total  intensity,  “/a”  and  “ Ib ”  are  the  inten¬ 
sities  of  the  voxels  containing  purely  tissue  “a”  or  tissue 
“b,”  respectively.  The  total  number  of  voxels  in  this 
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Table  1 

Spin  Echo  Volume  (in  Voxels)  and  the  Diameter  (in  mm)  Calculated  from  Spin  Echo  Volume  for  Each  Glass  Bead 


Bead 

1 

2 

3 

4 

5 

6 

7 

8 

9 

Measure 

Spin  echo  volume 

35.4 

37.5 

38.2 

39.1 

96.5 

103.8 

113.6 

516.6 

912.1 

Spin  echo  diameter 

2.0 

2.1 

2.1 

2.1 

2.9 

2.9 

3.0 

5.0 

6.0 

Actual  diameter 

2.0 

2.0 

2.0 

2.0 

3.0 

3.0 

3.0 

5.0 

6.0 

volume  of  interest  is  denoted  by  “IV,”  and  the  number  of 
voxels  occupied  by  tissue  “a”  is  denoted  by  “na.”  Conse¬ 
quently,  the  number  of  voxels  occupied  by  tissue  “b” 
can  be  expressed  as  N  —  na. 

By  varying  the  size  of  the  volume  of  interest,  the  total 
intensity  is  linearly  dependent  on  the  number  of  voxels 
in  the  volume  of  interest.  Although  “7b”  can  be  deter¬ 
mined  as  the  slope  in  the  fit  to  Eq.  7,  na  can  be  calcu¬ 
lated  from  the  intercept  if  “7a”  is  given  (na  may  not  be  an 
integer  as  partial  volume  is  included).  In  this  study,  “7b” 
corresponds  to  the  intensity  of  a  voxel  composed  purely 
of  gel,  whereas  “7a”  corresponds  to  the  intensity  of  a 
voxel  composed  purely  of  air  or  glass.  For  a  relatively 
large  air  bubble  or  glass  bead,  “7a”  is  dominated  by  the 
thermal  noise,  which  can  be  approximated  as  1.25  x 
crmag-gei»  where  amag_gei  is  the  measured  standard  devia¬ 
tion  of  the  gel  region  in  the  magnitude  images  (19).  For 
an  air  bubble  or  glass  bead  with  a  radius  generally  <3 
pixels,  “7a”  is  a  combination  of  thermal  noise  and  Gibbs 
ringing.  To  best  account  for  these  fluctuations,  “7a”  is 
calculated  from: 

nq  '  a mag— air  T  w2  "  1.25  ■  (Tmag— geL  for  air  bubble 
'  ^rnag  glass  T  w2  '  1*25  ■  rrmag_ger,  for  glass  bead 

[8] 

where  amag_ail.  and  amag_giass  are  the  measured  mean  values 
inside  the  bubble  and  glass  bead,  respectively;  w1  and  w2 
are  two  weighting  factors.  Based  on  our  simulations 
(explained  below),  wq  and  w2  were  empirically  deter¬ 
mined  from  simulations  to  be  0.4  and  0.6,  respectively,  to 
minimize  the  error  in  estimation  of  the  true  volume. 

Error  in  Volume  Measurement 

Although  a  regression  method  is  used  to  measure  the 
spin  echo  volume,  it  is  still  affected  by  partial  volume 
effects,  Gibbs  ringing,  as  well  as  random  noise.  The 
simulated  magnitude  images  at  TE  =  0  were  used  to 
mimic  spin  echo  magnitude  images  and  to  study  the 
error  in  spin  echo  volume  estimation.  In  addition,  to 
examine  the  stability  of  this  method  relative  to  thermal 
noise,  the  volume  measurement  evaluation  was  per¬ 
formed  10  times  for  each  simulated  sphere,  with  inde¬ 


pendently  generated  random  noise  for  each  of  these  sim¬ 
ulations.  The  errors  were  determined  by  comparing  the 
measured  volume  with  the  true  volume.  Note  that  this 
error  estimation  does  not  apply  for  the  apparent  volume 
which  is  determined  directly  from  the  binary  masks. 

RESULTS 

Simulations 

Magnetic  moments  for  simulated  spheres  were  calculated 
with  the  measured  susceptibilities  and  the  apparent  vol¬ 
ume  for  each  sphere  at  a  given  TE.  The  results  across  dif¬ 
ferent  TEs  are  shown  in  Fig.  1.  The  measured  volumes  at 
different  TEs  were  normalized  to  the  volume  at  the  lon¬ 
gest  TE,  whereas  the  measured  magnetic  moments  were 
normalized  to  the  true  magnetic  moment,  which  is  the 
product  of  input  volume  (i.e.,  the  true  volume)  of  the 
sphere  and  the  input  susceptibility  (true  susceptibility) 
9.4  ppm.  The  normalized  magnetic  moment  is  roughly  a 
constant  for  all  spheres.  However,  for  the  sphere  with  a 
radius  <2  pixels,  the  magnetic  moments  measured  in  the 
short  TE  range  have  more  fluctuations  than  those  meas¬ 
ured  at  longer  TEs.  In  addition,  the  magnetic  moments 
are  under-estimated  for  all  spheres.  The  mean  normal¬ 
ized  magnetic  moments  were  measured  as:  0.85  ±  0.04 
(radius  =  1  pixel),  0.82  ±  0.05  (radius  =  1.5  pixels),  0.81 
±  0.03  (radius  =  2  pixels),  and  0.81  ±  0.02  (radius  =  3 
pixels). 

After  the  magnetic  moments  were  obtained,  the  sus¬ 
ceptibility  values  were  corrected  using  the  actual  known 
volume  (i.e.,  true  volume)  using  Eq.  6.  Specifically,  the 
corrected  susceptibilities  are  7.95  ±  0.38  ppm  (radius  = 
1  pixel),  7.70  ±  0.43  ppm  (radius  =  1.5  pixels),  7.62  ± 
0.27  ppm  (radius  =  2  pixels),  and  7.63  ±  0.21  ppm 
(radius  =  3  pixels).  There  is  still  a  15-19%  underestima¬ 
tion  in  the  averaged  susceptibility  after  attempting  to 
correct  the  volume  of  the  sphere. 

To  evaluate  the  stability  of  the  volume  measuring 
method,  we  carried  out  10  simulations  for  each  sphere  at 
TE  =  0.  The  means  and  standard  deviations  of  the  per¬ 
centage  errors  relative  to  true  volume  for  each  sphere  are 
18.02  ±  27.26%  (radius  =  1  pixel),  1.89  ±  12.18%  (ra¬ 
dius  =  1.5  pixels),  3.67  ±  8.91%  (radius  =  2  pixels),  and 
2.09  ±  2.54%  (radius  =  3  pixels).  The  algorithm  failed  to 


Table  2 

Spin  Echo  Volume  (in  Voxel)  of  the  14  Air  Bubbles 


Bubble 

1 

2 

3 

4 

5 

6 

7 

Volume 

3.3 

15.1 

28.7 

42.2 

43.9 

82.8 

87.7 

Bubble 

8 

9 

10 

11 

12 

13 

14 

Volume 

92.7 

118.2 

170.5 

238.6 

288.7 

322.7 

897.2 
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FIG.  2.  Axial,  sagittal,  and  coronal  views  of  the  susceptibility  maps  with  TE  =  3.93  ms  (a,  b,  and  c)  and  TE  =  26.61  ms  (d,  e,  and  f). 
The  main  field  direction  is  in  “y"  direction.  Glass  bead  No.  9  in  Table  1  is  pointed  by  the  white  arrows.  The  air  bubbles  are  indicated  by 
the  white  dashed  arrows. 


quantify,  in  two  of  the  10  simulations  for  sphere  with  ra¬ 
dius  of  1  pixel.  Larger  errors  and  more  variations  of  the 
volume  measurements  were  seen  in  spheres  with  radii 
<2  pixels.  For  the  sphere  with  a  radius  of  3  pixels,  the 
error  in  the  volume  estimation  appears  to  be  within  5% 
using  the  proposed  method.  As  can  be  expected,  when 
the  object  radius  is  only  1  pixel,  the  volume  measure¬ 
ment  becomes  unstable. 

Phantom  Experiments 

A  total  of  14  air  bubbles  and  9  glass  beads  were  exam¬ 
ined  in  the  phantom  data.  The  measured  spin  echo  vol¬ 
umes  of  the  glass  beads  and  air  bubbles  are  summarized 
in  Tables  1  and  2.  In  these  two  tables,  the  glass  beads  as 
well  as  air  bubbles  are  sorted  based  on  their  spin  echo 
volumes,  from  small  to  large  objects.  The  diameters  of 
these  glass  beads  calculated  from  their  spin  echo  vol¬ 
umes  agree  reasonably  well  with  their  physically  mea¬ 
sured  diameters,  as  summarized  in  Table  1.  Also  note 
that  the  error  in  volume  measurement  is  unreliable  for 
spherical  objects  with  radii  <1.5  pixels  (14.13  voxels  for 
the  volume).  The  error  is  generally  larger  than  20%,  as 
shown  in  the  simulations.  Thus,  the  first  two  smallest 
air  bubbles  were  excluded  from  the  analysis. 

Figure  2  shows  three  orthogonal  views  of  the  suscepti¬ 
bility  map  of  the  largest  glass  bead  for  the  shortest  TE 


and  the  longest  TE.  Using  Eq.  6,  the  measured  suscepti¬ 
bilities  can  be  corrected  with  the  volume  estimated  from 
the  spin  echo  images.  These  results  are  summarized  in 
Table  3.  The  mean  of  the  corrected  susceptibility  values 
of  the  glass  beads  averaged  over  all  the  TEs  is  —1.82  ± 
0.17  ppm,  which  is  within  the  range  of  the  measured 
values  in  the  previous  study  (15).  The  mean  of  the  cor¬ 
rected  susceptibility  values  of  the  air  bubbles  is  6.66  ± 
0.85ppm.  This  is  to  be  compared  to  the  actual  suscepti¬ 
bility  of  9.4  ppm. 

DISCUSSION 

The  susceptibility  mapping  technique  using  the  regular¬ 
ized  Fourier-based  method  has  certain  advantages  over 
other  methods,  especially  in  terms  of  time  efficiency  and 
simplicity.  However,  it  suffers  from  problems  caused  by 
the  intrinsic  singularities  in  the  inverse  of  the  Green’s 
function,  as  well  as  partial  volume  effects  which  disrupt 
the  true  phase  behavior.  For  small  objects,  susceptibility 
quantification  using  the  inverse  method  (1)  yields  a  signif¬ 
icant  underestimation  of  the  susceptibility.  The  increased 
apparent  volume  at  long  TE  can  be  utilized  to  create  a 
larger  virtual  object  for  which  the  actual  susceptibility  can 
be  more  accurately  measured  and  thus  the  Fourier-based 
susceptibility  quantification  gives  a  relatively  smaller 
error  for  the  magnetic  moment.  At  this  point,  the 


Table  3 

Mean  Measured  and  Corrected  Susceptibilities  (in  ppm)  of  the  Glass  Beads  and  Air  Bubbles  at  Different  TEs 


TE/ms 

Glass  Bead 

Air  bubble 

Measured 

Corrected 

Measured 

Corrected 

3.93 

-1.50  ±  0.07 

-1.79  ±  0.13 

3.15  ±  1.16 

6.13  ±  0.77 

9.60 

-1.07  ±  0.32 

-1.76  ±  0.13 

1.70  ±  0.61 

6.44  ±  0.90 

15.27 

-0.86  ±  0.32 

-1.82  ±  0.19 

1.30  ±  0.45 

6.64  ±  0.76 

20.94 

-0.68  ±  0.28 

-1.82  ±  0.18 

1.07  ±  0.36 

6.88  ±  0.77 

26.61 

-0.59  ±  0.25 

-1.88  ±  0.21 

0.93  ±  0.34 

7.22  ±  0.74 
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susceptibility  close  to  the  actual  value  can  be  extracted 
from  the  estimated  magnetic  moment  with  an  estimation 
of  the  true  volume,  either  if  it  is  known  ahead  of  time,  or  it 
can  be  estimated  from  a  high-resolution  spin  echo  data  set. 

Based  on  the  discussions  above,  the  error  8Ay  in  the 
corrected  susceptibility  Ay  comes  from  the  estimated 
magnetic  moment  pa  =  Ax'F"  and  estimated  volume  (V). 
Through  error  propagation,  the  error  in  the  corrected  sus¬ 
ceptibility  is  given  by: 


As  we  can  see,  the  smaller  the  error  in  the  estimated  vol¬ 
ume,  the  smaller  the  error  in  the  corrected  susceptibility. 
This  equation  explains  the  error  seen  in  the  corrected 
susceptibility  of  the  air  bubbles  as  well  as  glass  beads. 

In  simulations,  where  the  true  volume  is  known,  the 
remnant  underestimation  in  the  averaged  corrected  sus¬ 
ceptibility  ranges  from  15  to  19%.  As  there  is  no  error  in 
the  true  volume,  this  error  must  be  owing  to  the  error  in 
the  apparent  volume  measurement  and  Ay'  quantification 
owing  to  the  regularization  process.  The  level  of  under¬ 
estimation  is  related  to  the  threshold  value  in  the  inverse 
of  the  Green’s  function.  A  smaller  threshold  leads  to  less 
underestimation,  but  more  streaking  artifacts  in  the  sus¬ 
ceptibility  maps.  The  regularized  Fourier-based  method, 
with  threshold  value  of  0.1,  can  lead  to  an  underestima¬ 
tion  of  around  13%  for  objects  with  radii  larger  than  3 
pixels  and  even  worse  for  smaller  objects  (1,6).  This  can 
be  viewed  as  a  systematic  error. 

In  phantom  studies,  after  using  the  volume  estimated 
from  the  spin  echo  data,  the  corrected  susceptibility  val¬ 
ues  of  the  air  bubbles  have  a  maximum  underestimation 
close  to  44%,  compared  to  the  theoretical  value  9.4  ppm. 
This  is  essentially  a  consequence  of  error  in  the  spin 
echo  volume  measurement  and  the  underestimation  of 
Ay'  quantified  using  the  Fourier-based  method.  To  over¬ 
come  these  limitations,  one  has  to  go  to  high-resolution 
images  that  can  minimize  volume  quantification  error 
and  to  relatively  longer  TEs  that  can  improve  accuracy 
in  the  magnetic  moment  quantification.  However,  the 
decreased  signal-to-noise  ratio  in  high-resolution  spin 
echo  images  may  introduce  additional  variation/noise  in 
the  final  volume  results. 

There  are  a  number  of  limitations  to  this  study. 
Although  a  forward  calculation  was  carried  out  to 
reduce  the  geometry  induced  field  variation,  remnant 
background  field  variation  still  exists.  To  best  account 
for  it,  the  phase  inside  the  spherical  objects  was  set  to 
the  local  average  phase.  This  also  helps  to  reduce  the 
large  variation  in  susceptibility  estimate  induced  by 
Gibbs  ringing  and  thermal  noise.  However,  this  phase 
correction  process  is  based  on  the  assumption  that  the 
object  of  interest  is  a  sphere.  For  nonspherical  objects, 
this  phase  correction  process  may  lead  to  variations  of 
magnetic  moment  at  different  TEs.  In  addition,  phase 
correction  also  creates  a  virtually  larger  object.  It  is  pos¬ 
sible  that  the  center  of  the  created  object  deviates  from 
the  true  center  of  the  original  object  of  interest.  This 
leads  to  additional  errors  even  for  spherical  objects,  as 


seen  from  simulations.  Thus,  the  phase  inside  the  spher¬ 
ical  object  has  significant  effects  to  this  method.  Theo¬ 
retically,  only  when  the  center  of  a  simulated  large 
sphere  coincides  with  the  original  center  of  the  sphere, 
and  when  the  background  phase  value  is  0,  we  can 
obtain  a  constant  magnetic  moment  across  different  TEs. 
Hence,  slight  variations  in  object  definition  from  a  bi¬ 
nary  mask,  which  is  used  for  phase  substitution,  can 
introduce  variations  in  magnetic  moment  values.  This  is 
the  essential  source  of  shape  dependence  of  the  pro¬ 
posed  method. 

Although  the  estimated  magnetic  moments  of  the  glass 
beads  are  almost  a  constant  over  different  TEs,  as  indi¬ 
cated  by  the  corrected  susceptibilities,  the  estimated 
magnetic  moments  of  the  air  bubbles  are  usually  larger  at 
a  longer  TE  than  at  a  shorter  TE.  This  can  be  understood 
by  the  fact  that  the  air  bubbles  are  not  perfect  spherical 
objects  compared  to  the  glass  beads.  In  fact,  most  of  the 
air  bubbles  have  ellipsoidal  shapes,  and  any  attempt  of 
phase  correction  inside  the  bubble  based  on  the  assump¬ 
tion  of  the  spherical  shape  will  cause  errors  in  the  sus¬ 
ceptibility  measurement  and  thus  lead  to  errors  in  the 
measurement  of  the  magnetic  moments. 

Generally  speaking,  for  small  objects  which  can  be 
well  approximated  as  spheres,  the  theoretically  expected 
errors  in  the  estimated  magnetic  moment  measurements 
are  within  20%  of  the  expected  values  and  can  be  fur¬ 
ther  reduced  by  adjusting  the  regularization  thresholds 
in  the  susceptibility  mapping  method.  Practically,  the 
errors  might  be  larger  owing  to  the  limited  knowledge  of 
the  true  volume.  Although  most  small  microbleeds  can 
be  well  approximated  as  spheres,  the  use  of  more  accu¬ 
rate  volume  estimation  methods  has  the  potential  to 
reduce  the  error  in  susceptibility  quantification  of 
microbleeds. 

CONCLUSIONS 

In  conclusion,  we  have  shown  that  for  very  small  struc¬ 
tures,  obtaining  accurate  magnetic  susceptibility  values 
is  limited  by  the  errors  in  the  volume  estimations  of 
these  structures  and  in  the  Fourier-hased  method  itself. 
Despite  this  inability  to  estimate  the  actual  volume  of  a 
small  object  accurately  (whether  it  is  an  air  bubble  or 
microbleed),  the  estimated  magnetic  moment  is  almost  a 
constant  over  different  TEs.  This  demonstrates  that  it  is 
possible  to  measure  the  magnetic  moment  at  a  longer  TE 
when  the  apparent  volume  is  increased  owing  to  T2* 
dephasing.  By  measuring  or  knowing  a  priori  the  actual 
volume  of  an  object,  it  is  possible  to  obtain  a  reasonable 
estimate  of  the  susceptibility. 
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Improving  Susceptibility  Mapping  Using  a 
Threshold-Based  K-Space/Image  Domain  Iterative 
Reconstruction  Approach 

J.  Tang,1  S.  Liu,1  J.  Neelavalli,2  Y.  C.  N.  Cheng,2  S.  Buch,1  and  E.  M.  Haacke1-3* 


To  improve  susceptibility  quantification,  a  threshold-based  k- 
space/image  domain  iterative  approach  that  uses  geometric 
information  from  the  susceptibility  map  itself  as  a  constraint  to 
overcome  the  ill-posed  nature  of  the  inverse  filter  is  intro¬ 
duced.  Simulations  were  used  to  study  the  accuracy  of  the 
method  and  its  robustness  in  the  presence  of  noise.  In  vivo 
data  were  processed  and  analyzed  using  this  method.  Both 
simulations  and  in  vivo  results  show  that  most  streaking  arti¬ 
facts  inside  the  susceptibility  map  caused  by  the  ill-defined 
inverse  filter  were  suppressed  by  the  iterative  approach.  In 
simulated  data,  the  bias  toward  lower  mean  susceptibility  val¬ 
ues  inside  vessels  has  been  shown  to  decrease  from  around 
10%  to  2%  when  choosing  an  appropriate  threshold  value  for 
the  proposed  iterative  method.  Typically,  three  iterations  are 
sufficient  for  this  approach  to  converge  and  this  process  takes 
less  than  30  s  to  process  a  512  x  512  x  256  dataset.  This 
iterative  method  improves  quantification  of  susceptibility  inside 
vessels  and  reduces  streaking  artifacts  throughout  the  brain 
for  data  collected  from  a  single-orientation  acquisition.  This 
approach  has  been  applied  to  vessels  alone  as  well  as  to  ves¬ 
sels  and  other  structures  with  lower  susceptibility  to  generate 
whole  brain  susceptibility  maps  with  significantly  reduced 
streaking  artifacts.  Magn  Reson  Med  69:1396-1407,  2013. 
©2012  Wiley  Periodicals,  Inc. 
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to  utilize  this  concept  in  MRI.  Unfortunately,  this 
inverse  process  is  ill-posed  and  requires  a  regularization 
procedure  to  estimate  the  SM.  There  are  a  variety  of 
approaches  to  tackle  this  problem  (7-24).  One  unique 
method  uses  a  multiple  orientation  data  acquisition  to 
remove  the  singularities  (17).  Constrained  regularizations 
(14,20,22,23)  have  shown  good  overall  results,  but  they 
require  longer  reconstruction  times  and  assumptions 
about  the  contrast  in  or  near  the  object  to  be  detected. 
Threshold-based,  single-orientation  regularization  meth¬ 
ods  (TBSO)  (11,15,18,24)  provide  the  least  acquisition 
time  and  the  shortest  computational  time  to  calculate 
SM.  However,  their  calculated  SMs  lead  to  underesti¬ 
mated  susceptibility  values  (Ay)  and  display  severe 
streaking  artifacts  especially  around  structures  with  large 
Ay,  such  as  veins  or  parts  of  the  basal  ganglia. 

Based  on  TBSO  approaches,  we  propose  an  iterative 
method  to  overcome  the  singularities  in  the  inverse  filter 
and  produce  improved  accuracy  for  susceptibility  map¬ 
ping.  In  this  approach,  we  iteratively  replace  k-space  val¬ 
ues  associated  with  the  SM,  y(k),  near  the  singularities 
to  obtain  an  almost  artifact  free  SM,  y(r).  The  k-space 
values  used  for  substitutions  are  estimated  using  struc¬ 
tural  information  from  the  masked  version  of  y(r).  Simu¬ 
lations  using  2D  cylinders  and  full  brain  3D  models 
were  performed  to  examine  the  efficacy  of  this  iterative 
approach.  High  resolution  human  data  are  also 
evaluated. 


Susceptibility  weighted  imaging  (SWI)  using  phase  infor¬ 
mation  has  become  an  important  clinical  tool  (1-3). 
However,  the  use  of  phase  information  itself  has  stimu¬ 
lated  great  interest  both  as  a  source  of  contrast  (4-6)  and 
a  source  for  producing  susceptibility  maps  (SM)  (7-24). 
The  impetus  for  solving  the  inverse  problem  from  mag¬ 
netic  field  perturbation  came  from  the  work  described  by 
Deville  et  al.  (25).  This  was  noted  by  Marques  and  Bow- 
tell  in  2005  (26).  Salomir  et  al.  (27)  were  the  first  group 


METHODS 

Briefly,  the  expression  for  the  susceptibility  distribution 
(26,27)  derived  from  the  phase  data  can  be  written  as 
(for  a  right  handed  system  (28)): 


y(r)  =  FT- 


1  x  FT 

$(r)  1 

[*(*) 

■ yB0TE\ 

where, 
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k  2 

g(k)  =- - ^ — %- 

3  kx2+k2 


[2] 


and  4>(r)  is  the  phase  distribution,  TE  is  the  echo  time,  -y 
is  the  gyromagnetic  ratio  for  hydrogen  protons,  B0  is  the 
main  field  strength,  kx,  ky,  and  kz  are  coordinates  in  k- 
space,  and  g[k)  is  the  Green’s  function  or  filter.  Clearly, 
the  analytic  inverse  filter  g_1(k)  =  l/g{k),  is  ill-posed 
when  g[k)  is  equal  or  close  to  zero,  i.e.,  points  on  or  near 
two  conical  surfaces  in  k-space  at  the  magic  angles  of 
54.7°  and  125.3°  from  the  B0  axis.  This  ill-posedness 
leads  to  severe  artifacts  (including  severe  streaking)  in 
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FIG.  1 .  Illustration  of  the  iterative  reconstruction  algorithm  to  obtain  artifact  free  x(r)  maps. 


x(k)  and  noise  amplification  (29).  Thus,  for  a  proper  pixel- 
by-pixel  reconstruction  of  xW,  recovering  the  correct  val¬ 
ues  of  xW)  near  the  region  of  singularities  is  critical. 


K-space  Iterative  Approach 

If  the  shapes  of  the  structures  of  interest  are  known,  then 
one  can  use  this  information  in  the  SM  to  create  a  more 
accurate  k-space  of  said  SM  in  the  conical  region.  The 
structure  of  the  vessels  is  obtained  directly  from  the  first 
pass  SM  x/— oM-  The  detailed  steps  of  the  iterative 
method  are  discussed  below  and  shown  in  Fig.  1. 

•  Step-1:  An  initial  estimate  of  the  SM,  X/=oW>  is 
obtained  by  applying  a  regularized  version  of  the 
threshold-based  inverse  filter,  gfeg(k)  (18),  in  Eq.  1 
using  the  suggested  threshold  value,  thr  =  0.1.  The 
subscript  “i”  denotes  the  SM  after  the  ith  iteration 
(“i  =  0”  denotes  the  initial  step  before  doing  the 
iterative  method  and  i  =  1  for  the  first  iteration). 

•  Step-2:  The  geometry  of  the  structures  of  interest  is 
extracted  from  Xj=oW  using  a  binary  vessel  mask, 
i.e.  outside  the  veins,  the  signal  in  the  mask  is  set  to 
zero,  and  inside  it  is  set  to  unity.  Since  streaking 
artifacts  associated  with  veins  in  the  SM  are  usually 
outside  the  vessels,  after  multiplying  the  X;M  map 
by  the  mask,  little  streaking  remains  in  the  SM.  This 
leads  to  Xvm.iW  as  shown  in  Fig.  lb.  Vessel  mask 
generation  will  be  addressed  in  the  next  section. 

•  Step-3:  Xvm.iW)  is  obtained  by  Fourier  transforma¬ 
tion  Of  Xvm.iW  (Fig.  lc). 

•  Step-4:  The  predefined  ill-posed  region  of  k-space 
in  Xvm.iW)  is  extracted  (Fig.  Id).  These  extracted  k- 
space  data  are  denoted  by  Xvm.cone.iW)-  The  size  of 
Xvm.cone.iW)  is  decided  by  a  threshold  value,  a, 
which  is  assigned  to  g(k).  For  the  matrix  size  512  x 
512  x  512,  the  percentages  of  the  cone  region  in 


k-space  for  a  given  a,  are  2.4%  (a  =  0.01),  24.1%  (a 
=  0.1),  47.1%  (a  =  0.2),  and  70.6%  (a  =  0.3),  respec¬ 
tively.  When  a  increases,  the  size  of  Xvm.cone.iW) 
increases  too.  If  a  is  increased  too  much  then  most 
of  the  original  information  will  be  lost. 

•  Step-5:  Data  from  Xvm.cone.iW)  and  Xi=oW)  (Fig.  le) 
are  merged.  This  means  part  of  Xi=oW)  has  been 
replaced  by  Xvm.cone.iW)-  The  merged  data  are 
denoted  by  Xmerged,iW)  (Fig.  If). 

•  Step-6:  Inverse  Fourier  transformation  of  XmergediW) 
gives  the  improved  SM,  Xi+iW  (Fig.  lg). 

•  Step-7:  XiW  in  step-1  is  replaced  by  Xi+iW  from 
step-6  and  the  algorithm  is  repeated  until 

[(x;-W-xmW)2]/iv  <  e  [3] 

where  N  is  the  number  of  pixels  in  x;W  and  8  is  the 
tolerance  value  (chosen  here  to  be  0.004  ppm). 

Binary  Vessel  Mask  Generation 

The  binary  vessel  mask  was  generated  using  threshold¬ 
ing  from  the  xW  map  itself.  The  detailed  steps  are  dis¬ 
cussed  below  and  shown  in  Fig.  2. 

•  Step-1:  A  threshold,  tha,  is  applied  to  Xr=oW  to  cre¬ 
ate  an  initial  binary  vessel  mask,  M0.  The  pixels 
whose  susceptibility  values  are  lower  than  thi  will 
be  set  to  zero  while  those  greater  than  or  equal  to 
tha  will  be  set  to  unity.  In  this  study,  a  relatively 
low  susceptibility  of  0.07  ppm  is  used  for  tht  to  cap¬ 
ture  most  vessels.  However,  this  choice  of  threshold 
inevitably  includes  other  brain  structures  in  M0  that 
have  high  susceptibility. 

•  Step-2:  A  morphological  closing  operation  is  per¬ 
formed  to  fill  in  holes  in  M0  to  generate  an  updated 
mask  Mi. 
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FIG.  2.  Illustration  of  the  binary  vessel  mask  generation  process. 


•  Step-3:  A  median  filter  is  applied  to  remove  noise 
in  Ma  and  create  M2. 

•  Step-4:  False  positive  data  points  from  M2  are 
removed  as  follows:  First,  the  xM  map  is  Mipped 
over  five  slices  centered  about  the  slice  of  interest  to 
better  obtain  contiguous  vessel  information,  as  seen 
in  XmipW-  Second,  another  threshold,  th2  =  0.25 
ppm,  is  performed  on  XmipM  to  create  a  new 
Xmip_viiiM  and  binary  mask  MP,  which  only  contains 
predominantly  vessels.  Here,  0.25  ppm  was  chosen 
to  isolate  the  major  vessels  in  the  MIP  image.  Third, 
each  slice  from  M2  is  compared  with  MP  on  a  pixel- 
by-pixel  basis  to  create  M3.  If  a  data  point  from  M2 
does  not  appear  on  MP,  this  data  point  will  be 
treated  as  a  false  positive  and  removed  from  M2,  oth¬ 
erwise  this  point  is  retained.  This  process  can  be 
equally  well  applied  to  extract  other  tissues  by 
choosing  appropriate  values  for  tha  and  th2. 

2D  Cylinder  Simulations 

Simulation  of  a  two  dimensional  cylinder  and  its 
induced  phase  was  first  performed  using  a  8192  x  8192 
matrix.  A  lower  resolution  complex  image  was  then 
obtained  by  taking  the  Fourier  transform  of  this  matrix 
and  applying  an  inverse  Fourier  transform  of  the  central 
512  x  512  matrix  in  k-space.  This  procedure  is  to  simu¬ 
late  Gibbs  ringing  effects  caused  by  finite  sampling 
which  we  usually  see  in  conventionally  required  MR 
data  sets.  Gibbs  ringing  comes  from  discontinuities  in 
both  the  magnitude  and  phase  images.  To  avoid  Gibbs 
ringing  from  magnitude  discontinuities,  we  used  a  mag¬ 
nitude  image  with  a  uniform  signal  of  unity.  Cylinders 
with  diameters  32,  64,  128,  256,  512,  and  1024  were 
simulated  on  8192  x  8192  matrices  and  their  effective 


diameters  were  2,  4,  8,  16,  32,  and  64  on  512  x  512  mat¬ 
rices.  All  phase  simulations  were  performed  using  a  for¬ 
ward  method  (8,26,27,30)  with  B0  =  3  T,  Ay  =  0.45  ppm 
in  SI  units,  TE  =  5  ms,  and  the  cylinder  perpendicular  to 
the  main  magnetic  field.  The  susceptibility  value  of  0.45 
ppm  represents  venous  blood  when  the  hematocrit  (Hct) 
=  0.44,  Axdo  =  4tt  0.27  ppm  (31)  and  the  oxygen  satura¬ 
tion  level  =  70%,  where  Axdo  is  the  susceptibility  differ¬ 
ence  between  fully  deoxygenated  and  fully  oxygenated 
blood  (32).  A  relatively  short  echo  time  was  chosen  to 
avoid  phase  aliasing  that  can  affect  the  estimated  suscep¬ 
tibility  values. 

Selection  of  a  TBSO  Method  to  Generate  the  Xi=oM  Map 

TBSO  methods  (11,15,18,24)  use  a  truncated  g[k)  to  solve 
the  singularity  problem  in  the  inverse  filter  g_1(k)  when 
g(k)  is  less  than  a  predetermined  threshold  value,  thr. 
When  g{k)  <  thr,  g_1(k)  is  either  set  to  zero  (11,24);  or  to 
1/thr  (15);  or  set  to  g_1(k)  =  1/thr  first  and  then  g^1(k)  is 
brought  smoothly  to  zero  as  k  approaches  kzo.  This 
smoothing  is  accomplished  by  multiplying  g_1(k)  by 
a2(kz)  with  a(kz)  =  (kz  -  kzo)/  I  kzthr  -  kzo  I  where  kz  is 
the  z  component  of  that  particular  point  in  k-space,  kzo 
is  the  point  at  which  the  function  g_1(k)  becomes  unde¬ 
fined,  and  kzthr  is  the  kz  coordinate  value  where  I  g(k)  I 
=  thr  (18). 

SMs  using  the  methods  in  Refs.  11,15,18  were  calcu¬ 
lated  based  on  Eq.  1  using  the  2D  cylindrical  model. 
Equation  1  can  be  used  to  calculate  the  SM  for  the  simu¬ 
lated  2D  cylinder  model  perpendicular  to  the  main 
field  since  the  2D  perpendicular  model  is  a  special  case 
of  the  3D  model  with  1  slice  [10].  Streaking  artifacts  are 
obvious  in  all  three  SMs  (figures  are  not  shown).  The 
calculated  mean  susceptibility  values  inside  the  cylinder 
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FIG.  3.  a:  A  transverse  view  of  the  3D  brain  model,  b:  The  simulated  phase  map  from  the  model  using  parameters:  B0  =  3  T  and  TE  = 
18  ms  which  are  consistent  with  imaging  parameters  in  the  real  data  (c).  Images  (b)  and  (c)  have  the  same  window  level  setting. 


are  around  0.40  ±  O.Olppm  for  all  SMs.  The  background 
noise  levels,  (i.e.,  standard  deviation  of  the  susceptibility 
values)  measured  from  a  region  outside  the  streaking  ar¬ 
tifact  in  SM  using  Ref.  18  are  around  1/2  to  2/3  of  the 
background  noise  levels  in  SMs  using  Refs.  11,15  with 
thr  =  0.06,  0.07,  and  0.1,  which  are  the  optimal  thresh¬ 
old  values  suggested  in  Refs.  11,12,18.  Given  this  result, 
the  method  in  (18)  was  chosen  to  generate  a  Xi=oM  map. 

Finding  an  Optimal  Threshold  Value 

To  find  the  optimal  threshold,  a  series  of  xM  maps  were 
reconstructed  by  the  iterative  method  using  threshold 
values  a  of  0.01,  0.03,  0.07,  0.1,  0.15,  0.2,  0.25,  and  0.3. 
The  larger  this  threshold,  the  closer  the  final  estimate  for 
x(r)  will  be  to  XvmW-  The  optimal  threshold  value  was 
found  by  comparing  the  accuracy  of  the  estimated  sus¬ 
ceptibility  values  as  well  as  the  effects  on  reducing 
streaking  artifacts  in  the  reconstructed  xM  maps.  To 
study  the  effect  of  noise  in  xM  maps  due  to  the  noise  in 
phase  images,  complex  datasets  for  cylinders  of  diameter 
2,  4,  8,  32  voxels,  respectively,  were  simulated  with 
Gaussian  noise  added  to  both  real  and  imaginary  chan¬ 
nels.  Noise  was  added  in  the  complex  images  to  simulate 
a  SNRmagnitude  of  40:1,  20:1,  10:1,  and  5:1  in  the  magni¬ 
tude  images.  Since  o-p^ase  =  l/SNRmagnitude,  this  corre¬ 
sponds  to  crphase  =  0.025,  0.05,  0.1,  and  0.2  radian. 

To  estimate  the  improvement  in  the  SM  by  the  itera¬ 
tive  method,  we  used  a  root  mean  squared  error  (RMSE) 
to  measure  streaking  artifacts  outside  the  cylinder.  Back¬ 
ground  noise  in  the  SM  is  measured  in  a  region  away 
from  all  major  sources  of  streaking  artifacts  to  compare 
the  noise  measured  in  the  phase  image  (i.e.,  so  we  can 
correlate  noise  in  the  phase  with  the  expected  noise 
enhancement  from  the  inversion  process). 

Effect  of  High-Pass  Filter 

The  effect  of  high-pass  (HP)  filtering  the  phase  data  on 
the  xM  map  generated  by  the  iterative  method  was  also 
studied.  Phase  images  of  a  cylindrical  geometry  with 
diameters  of  2,  4,  8,  16,  32,  and  64  voxels  were  simu¬ 
lated.  Homodyne  HP  filters  (33)  with  a  2D  hanning  filter 
(full  width  at  half-maximum,  FWHM  =  4,  8,  16,  and  32 
pixels)  were  applied  on  these  phase  images  in  both  in¬ 


plane  directions.  SM  reconstructions  were  stopped  based 
on  the  criteria  in  step  7  of  the  iterative  process. 

Three  Dimensional  Brain  Model  Simulations 

To  address  the  potential  of  the  iterative  technique  to 
improve  the  SM  of  general  structures  such  as  the  basal 
ganglia,  a  3D  model  of  the  brain  was  created  including 
the:  red  nucleus  (RN),  substantia  nigra  (SN),  crus  cerebri 
(CC),  thalamus  (TH),  caudate  nucleus  (CN),  putamen 
(PUT),  globus  pallidus  (GP),  gray  matter  (GM),  white 
matter  (WM),  cerebrospinal  fluid  (CSF),  and  the  major 
vessels  (34).  The  structures  in  the  3D  brain  model  were 
extracted  from  two  human  3D  TV  weighted  and  T2 
weighted  data  sets.  Basal  ganglia  and  vessels  are  from 
one  person;  GM  and  WM  are  from  the  other  person’s 
data  set.  Since  all  structures  are  from  in  vivo  human 
data  sets,  this  brain  model  represents  realistic  shapes 
and  positions  of  the  structures  in  the  brain.  Susceptibil¬ 
ity  values  in  parts  per  million  (ppm)  for  the  structures 
SN,  RN,  PUT,  and  GP,  were  taken  from  Ref.  12  and 
others  were  from  measuring  the  mean  susceptibility 
value  in  a  particular  region  from  SMs  using  Ref.  18  from 
in  vivo  human  data:  RN  =  0.13,  SN  =  0.16,  CC  =  —0.03, 
TH  =  0.01,  CN  =  0.06,  PUT  =  0.09,  GP  =  0.18,  vessels  = 
0.45,  GM  =  0.02,  CSF  =  -0.014,  and  WM=0.  All  struc¬ 
tures  were  set  inside  a  512  x  512  x  256  matrix  of  zeros. 
The  phase  of  the  3D  brain  model  was  created  by  apply¬ 
ing  the  forward  method  (8,26,27,30)  to  the  3D  brain 
model  with  different  susceptibility  distributions  using 
the  imaging  parameters:  TE  =  5  ms  and  B0  =  3  T.  A  com¬ 
parison  between  the  phase  maps  from  this  brain  model 
and  a  real  data  set  is  shown  in  Fig.  3.  To  match  the 
imaging  parameters  of  the  real  data  set,  B0  =  3  T  and  TE 
=  18  ms  were  applied  for  the  results  presented  in  Fig.  3. 
Except  for  Fig.  3,  all  other  figures  in  the  paper  associated 
with  the  3D  brain  were  simulated  by  using  TE  =  5  ms. 

In  Vivo  MR  Data  Collection  and  Processing 

A  standard  high-resolution  3D  gradient  echo  SWI 
sequence  was  used  for  data  acquisition.  A  transverse  0.5 
mm  isotropic  resolution  brain  dataset  was  collected  at  3 
T  from  a  23-year-old  healthy  volunteer.  The  sequence  pa¬ 
rameters  were:  TR  =  26  ms,  flip  angle  =  15°,  read 
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FIG.  4.  Simulations  showing  the  comparison  of  the  calculated  susceptibility  distributions  for  a  cylinder  perpendicular  to  B0  at  different 
threshold  values  (a)  applied  to  g(k)  as  well  as  the  initial  x/  =  oW  map.  The  direction  of  B0  is  indicated  by  a  black  long  arrow.  The  suscep¬ 
tibility,  Ax.  inside  the  cylinder  is  0.45  ppm.  a:  The  comparison  of  the  converged  x/=o(r)  maP  with  the  x/=oM  map  for  a  diameter  of  32- 
pixel  cylinder,  where  b  is  the  iterative  step  required  to  reach  convergence.  In  this  data,  b  =  2  when  a  =  0.03,  b  =  3  when  a  =  0.1  and  b 
=  4  when  a  =  0.2  when  <rPhase  =  0.  The  top  row  of  images  shows  simulations  with  no  phase  noise.  The  second  and  the  third  row  show 
simulations  with  added  phase  noises  <jphase  =  0.025  and  0.05  radian,  respectively.  The  first  column  of  images  show  initial  xj-o(r)  maps 
for  reference,  b:  The  variation  of  the  mean  calculated  susceptibility  inside  the  cylinder  with  different  threshold  value,  a,  for  diameter  (d) 
=  2,  4,  8,  and  32  pixels  cylinders.  The  mean  susceptibility  value  is  independent  of  the  noise  level;  therefore,  only  mean  values  from 
'’■phase  =  0  were  provided,  c:  The  variation  of  the  RMSE  of  the  susceptibility  values  outside  the  cylinder  as  a  function  of  the  threshold 
value,  a,  and  the  noise  level.  The  d  =  32  pixels  cylinder  was  used  to  generate  (c).  The  range  of  the  gray-scale  bars  is  chosen  to  high¬ 
light  the  artifacts  in  the  images.  It  does  not  reflect  the  quantified  higher  susceptibility  values  inside  cylinders. 


bandwidth  =  121  Hz/pixel,  TE  =  14.3  ms,  192  slices,  and  a 
matrix  size  of  512  x  368.  To  reconstruct  X;=oM  with  mini¬ 
mal  artifacts,  the  following  steps  were  carried  out: 

1.  The  unwanted  background  phase  variations  were 
removed  using  either:  (a)  a  homodyne  HP  filter 
(FWHM  =  16  pixels)  (33)  or  (b)  Prelude  in  FMRIB 
Software  Library  (FSL)  (35)  to  unwrap  the  phase,  fol¬ 
lowed  by  the  process  of  Sophisticated  Harmonic  Arti¬ 
fact  Reduction  for  Phase  data  (SHARP)  (36)  with  a 
filter  radius  of  6  pixels.  To  reduce  artifacts  in  the  cal¬ 
culated  SMs,  regions  with  the  highest  phase  devia¬ 
tions  due  to  air/tissue  interfaces  were  removed 
manually  from  the  HP  filtered  phase  images  and  the 
phase  in  those  regions  were  set  to  zero. 

2.  A  complex  threshold  approach  (37)  was  used  to 
separate  the  brain  from  the  skull. 

3.  The  phase  image  with  an  original  matrix  size  of  512 
x  368  x  192  was  zero  filled  to  512  x  512  x  256  to 


increase  the  field-of-view  and  to  avoid  streaking 
artifacts  caused  by  the  edge  of  brain  to  alias  back  to 
the  reconstructed  SM. 

4.  The  regularized  inverse  filter,  g^glk)  (18)  was 
applied  to  obtain  x^oM,  followed  by  the  iterative 
process  using  a  =  0.1.  For  in  vivo  data,  the  iterative 
program  was  terminated  at  the  third  iterative  step. 

RESULTS 

Selection  of  Threshold  Level  Based  on  Simulations 

To  find  the  optimal  threshold  value,  SMs  were  recon¬ 
structed  using  a  =  0.01,  0.03,  0.07,  0.1,  0.15,  0.2,  0.25, 
and  0.3,  respectively,  with  different  noise  levels  (Fig.  4). 
The  streaking  artifacts  shown  in  X;=oM  (the  first  column 
in  Fig.  4a)  have  been  significantly  reduced  by  the  itera¬ 
tive  method  and  fall  below  the  noise  level  when  a  >= 
0.1.  Also,  when  a  >=  0.1,  the  mean  susceptibility  value 
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inside  the  cylinder  was  found  to  increase  to  0.44  ppm 
when  the  diameter  of  the  cylinder  was  larger  than  8  pix¬ 
els  (Fig.  4b)  and  this  trend  is  independent  of  the  object 
size  and  the  noise  in  the  phase  image.  The  optimal  result 
in  terms  of  obtaining  the  true  susceptibility  value  was 
with  a  threshold  of  0.1.  Figure  4c  shows  a  plot  of  RMSE 
of  the  susceptibility  values  from  the  whole  region  out¬ 
side  the  32-pixel  cylinder  using  different  a.  The  RMSE 
of  the  susceptibility  values  decreases  as  a  increases. 
Therefore,  for  vessels,  a  value  of  a  =  0.3  would  be  the 
optimal  value.  Flowever,  a  large  threshold  value  means 
replacing  more  original  k-space  with  the  k-space  only 
consisting  of  vessel  information  which  will  reduce  the 
signals  from  other  brain  structures  and  blur  these  struc¬ 
tures.  Since  the  SM  using  a  =  0.1  already  reveals  the 
optimal  susceptibility  value  for  the  vessels  and  an  ac¬ 
ceptable  RMSE,  it  is  appropriate  to  choose  0.1  for  more 
general  applications  to  study  the  entire  brain. 

Figure  4a  compares  the  converged  X/=z>(r)  map  with  the 
Xi=oM  map,  where  b  is  the  iterative  step  required  to 
reach  convergence.  In  this  data,  b  =  2  when  a  =  0.01 
and  0.03,  b  =  3  when  a  =  0.07,  0.1,  and  0.15  and  b  =  4 
when  a  =  0.2,  0.25,  and  0.3  for  Up^ase  =  0-  When  Uphase 
increases,  more  iterative  steps  were  required  to  reach 
convergence.  For  instance,  the  maximum  iterative  step 
number  is  9  when  uphase  =  0-2  radians.  Using  a  noise 
level  of  0.025  radian  in  the  phase  image  as  an  example, 
greg(k)  (18)  leads  to  a  susceptibility  noise  of  roughly 
0.025  ppm  in  the  Xj=oM  map.  The  iterative  approach 
leads  to  a  slight  decrease  in  background  noise,  0.021 
ppm,  in  X;=3(r)  map  when  a  =  0.1.  The  background  noise 
was  measured  in  a  region  outside  the  streaking  artifact 
indicated  by  the  black  circle  in  Fig.  4a.  The  overall 
decrease  in  RMSE  in  the  background  (Fig.  4c)  is  a  conse¬ 
quence  of  both  a  decrease  in  streaking  artifacts  and  a 
reduction  in  thermal  noise  contribution. 

Selection  of  the  Optimal  Iterative  Step 

The  inverse  process  (18)  was  applied  to  the  dipole  field 
in  Fig.  5a  to  give  the  X;=oM  map  shown  in  Fig.  5b; 
prominent  streaking  artifacts  are  evident  in  this  image. 
Streaking  artifacts  are  significantly  reduced  at  each  step 
of  the  iterative  method  quickly  reaching  convergence 
(Fig.  5c-e).  The  largest  improvement  is  seen  in  the  first 
iterative  step,  which  is  verified  by  Fig.  5f,  showing  the 
difference  between  Fig.  5c  (x;=iM  map)  and  Fig.  5b 
(Xi=oM  map).  After  the  second  iteration,  we  can  see 
some  minor  streaking  reductions  (Fig.  5g,  the  difference 
between  the  Xi=iM  map  and  Xi-'zM  map).  The  mean  sus¬ 
ceptibility  value  approaches  0.44  ppm  in  a  single  step. 
Similar  results  (not  shown)  are  also  obtained  when  the 
iterative  method  is  run  with  different  aspect  ratios 
between  the  in-plane  resolution  and  the  through  plane 
resolution  (such  as  1:2  and  1:4).  The  iterative  results 
always  lead  to  higher  final  susceptibility  values  com¬ 
pared  to  the  initial  value  in  Xi=oM-  Finally,  even  when 
an  HP  filter  is  applied,  up  to  a  10%  increase  in  the  sus¬ 
ceptibility  is  realized  (Fig.  5i).  The  SMs  of  large  vessels 
benefit  from  a  low  order  HP  filter  (FWMH  =  4  pixels) 
and  small  vessels  up  to  8  pixels  benefit  from  a  HP  filter 
(FWMH  =  16  pixels). 


Effect  of  the  Iterative  Approach  on  Surrounding  Brain 
Tissues  in  the  3D  Brain  Model 

SM  Reconstruction  Using  a  Vessel  Mask  Only 

Figure  6a, d  represents  X;=oM>  without  noise  and  with 
0.025  radians  of  noise  in  phase  images.  Figure  6f  is  the 
vessel  map.  Streaking  artifacts  (delineated  by  the  black 
arrows)  are  obvious  in  Fig.  6a, d  and  significantly 
reduced  in  the  Xr=3W  maps  (Fig.  6b, e)  using  a  =  0.1.  Fig¬ 
ure  6c  is  the  Xr=3(r)  map  using  a  =  0.2.  As  can  be  seen, 
when  a  increases,  the  iterative  method  still  works  for 
vessels,  but  brain  tissues  become  more  blurred.  Figure  7a 
plots  the  mean  susceptibility  values  inside  the  vessel 
(vein  of  Galen),  GP,  SN,  RN,  PUT,  and  CN  from  X;=3(r) 
maps  generated  by  using  a  =  0.1,  0.15,  0.2,  0.25,  and  0.3, 
respectively.  The  susceptibility  value  in  the  brain  model 
and  Xi=oM  map  are  also  provided  in  the  plot  as  referen¬ 
ces.  Generally,  the  susceptibility  values  of  brain  tissues 
except  vessels  decrease  as  a  increases  while,  for  vessels, 
the  susceptibility  value  is  0.41  ppm  in  the  Xj=oM  map 
and  is  increased  to  0.45  ppm  in  the  Xi=3(r)  maps. 

SM  Reconstruction  Using  a  Mask  Including  Vessels  and 
Brain  Structures 

The  iterative  method  is  not  limited  to  improving  SM 
from  just  vessels;  it  can  also  be  applied  to  the  entire 
brain.  Figure  6g  shows  a  coronal  view  of  the  Xi=oM  map 
for  the  brain  model.  The  Xi=3W  map  using  a  mask  keep¬ 
ing  all  major  structures  (GP,  SN,  RN,  PUT,  CN)  and  ves¬ 
sels  is  shown  in  Fig.  6h.  In  practice,  this  is  equivalent  to 
setting  thresholds  in  the  Xj=oM  map  to  be  greater  than 
0.09  ppm  to  extract  all  these  high  susceptibility  struc¬ 
tures  from  the  Xi=oM  map  to  create  the  mask.  Figure  6h 
reveals  that  streaking  artifacts  associated  with  veins  as 
well  as  all  major  structures  have  been  reduced.  Figure  6i 
shows  the  difference  between  Fig.  6g,  h.  In  addition, 
streaking  artifacts  sometimes  cause  the  appearance  of 
“false”  structures.  For  instance,  there  is  no  internal  cap¬ 
sule  (IC)  included  in  the  model  (Fig.  61),  yet  we  see  an 
IC  like  structure  in  the  X;=oM  map  (Fig.  6j)  (indicated  by 
a  dashed  white  arrow  in  Fig.  6j).  The  iterative  method 
removes  the  streaking  artifacts  and  the  “false”  IC  (Fig. 
6k).  Figure  7b  shows  susceptibility  values  in  each  struc¬ 
ture  in  the  brain  model  for  Xi=oM  and  Xi=3(r)  when  the 
mask  includes  vessels  and  all  major  structures.  The 
underestimated  susceptibility  values  of  all  major  struc¬ 
tures  and  vessels  in  the  X;=oM  map  have  been  recovered 
by  the  iterative  method  in  the  X;=3W  map. 

Effect  of  Errors  in  the  Vessel  Map 

Accurately  extracting  vessels  from  \i= oM  is  critical  for 
the  iterative  method.  Figure  8b-d  and  the  corresponding 
enlarged  views  (Fig.  8f-h)  from  the  rectangular  region 
indicated  in  Fig.  8a  show  the  Xr=3(r)  maps  using  an  accu¬ 
rate  (Fig.  8j),  a  dilated  (Fig.  8k),  and  an  eroded  (Fig.  81) 
vessel  map  to  show  the  effect  of  errors  in  the  vessel 
mask  on  the  \i=  sM  map.  The  dilated  and  eroded  vessel 
maps  were  generated  using  Matlab  functions  based  on  a 
3-by-3  square  structuring  element  object.  The  susceptibil¬ 
ity  values  measured  from  a  vein  indicated  by  the 
white  arrow  in  Fig.  8e  are  0.40  ±  0.03  ppm  (Fig.  8e), 
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FIG.  5.  a:  Phase  images  from  a  cylinder  with  a  diameter  of  32  pixels  are  simulated  with:  Ax  =  0.45  ppm,  S0  =  3  T,  and  TE  =  5  ms.  The 
cylinder  is  perpendicular  to  the  main  field.  No  thermal  noise  was  added  in  these  images,  b:  The  initial  xi-o(r)  map.  c:  The  SM  from  the 
first  iteration,  x/=iM  map,  (d)  x/=zW  map  and  (e)  x/=3(r)  maP  using  threshold  value  a  =  0.1.  The  SM  has  converged  at  x/= al/)  map.  The 
streaking  artifacts  are  reduced  as  the  number  of  iterative  steps  increases,  f:  The  difference  image  of  (c)  subtracted  from  (b)  illustrates 
that  the  streaking  artifacts  were  reduced  by  the  iterative  procedure  and  the  largest  improvement  happens  in  this  first  iterative  step,  g: 
The  difference  image  of  the  x/=iM  map  subtracted  from  the  x/=2 M  map  indicates  that  the  streaking  artifacts  were  further  reduced  by  the 
second  iterative  step,  h:  The  difference  image  of  x*  2W  map  subtracted  from  x/= al/)  map  shows  much  less  improvement  at  the  third  iter¬ 
ative  step.  Thus  it  indicates  a  convergence  of  the  iterative  procedure.  All  images  were  set  to  the  same  window  level  setting  for  direct 
comparisons  and  for  enhancing  the  presence  of  the  streaking  and  the  remnant  error,  i:  The  effect  of  the  iterative  approach  on  the 
changes  in  susceptibility  values  from  HP  filtered  phase  images.  Differences  between  the  values  in  iterative  and  non-iterative  susceptibil¬ 
ity  map  reconstruction  (i.e.,  XconvergedW  -  x;-oM)  from  HP  filtered  phase  images  are  plotted  for  different  filter  sizes.  Results  for  four  filter 
sizes  (FWHM  =  4,  8,  16,  and  32  pixels)  are  shown  here.  Applying  an  HP  filter  leads  to  an  underestimation  of  Ax  (18).  The  iterative 
approach  helps  to  improve  the  accuracy  of  the  estimated  susceptibility  values.  The  range  of  the  gray-scale  bars  is  chosen  to  highlight 
the  artifacts  in  the  images.  It  does  not  reflect  the  quantified  higher  susceptibility  values  inside  cylinders. 


0.45  ±  0.03  ppm  (Fig.  8f),  0.45  ±  0.03  ppm  (Fig.  8g),  and 
0.40  ±  0.07  ppm  (Fig.  8h),  respectively.  The  iterative 
method  still  works  if  the  vessel  is  slightly  enlarged  but 
does  little  to  change  the  original  Xi=oM  map  if  the  ves¬ 
sels  are  too  small  or  absent  in  the  mask.  As  we  just  dis¬ 
cussed,  streaking  artifacts  produced  “false”  vessels  indi¬ 
cated  by  the  dashed  black  arrow  in  Fig.  8e  since  these 
vessels  are  not  in  the  model  (Fig.  8i).  These  false  vessels 
disappeared  in  Fig.  8f. 

Results  from  the  In  Vivo  Dataset 

In  the  in  vivo  example,  we  compare  the  differences 
between  SHARP  (Fig.  9a-d)  and  a  homodyne  HP  filter 
(FWHM  =  16  pixels)  (Fig.  9e-h).  Compared  to  the  trans¬ 
verse  view,  streaking  artifacts  are  more  obvious  in  the 
sagittal  or  coronal  view.  Figure  9a  shows  the  Xi=oM  map 
with  severe  streaking  artifacts.  The  streaking  artifacts 


were  significantly  reduced  in  the  Xi=3W  map  (Fig-  9b) 
using  a  =  0.1.  The  streaking  artifacts  associated  with  the 
superior  sagittal  sinus  vein  (indicated  by  two  black 
arrows  in  Fig.  9a)  were  significantly  decreased  in  Fig. 
9b, d.  The  subtracted  image  (Fig.  9c),  Fig.  9a  minus  Fig. 
9b,  reveals  the  removed  streaking  artifacts.  These  streak¬ 
ing  artifacts  are  one  of  the  reasons  why  the  Xi=oM  maps 
appear  noisy.  In  the  Xj=3M  map,  the  reduction  in  streak¬ 
ing  artifacts  from  individual  veins  leads  to  a  decrease  of 
noise  therefore  an  increased  SNR  of  veins.  If  veins  are 
the  only  interest,  even  a  threshold  of  0.2  can  work  rea¬ 
sonably  well  (Fig.  9d).  Two  relatively  big  veins,  Vl  and 
V2,  indicated  by  the  white  dashed  and  white  solid 
arrows,  respectively,  in  Fig.  9b,  were  chosen  to  measure 
the  susceptibility  values.  Results  are  provided  in  Table 
1.  The  susceptibility  values  of  these  two  veins  have  been 
improved  by  roughly  16%  by  the  iterative  method.  The 
standard  deviation  of  the  susceptibility  values  measured 
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FIG.  6.  Results  before  and  after  the  iterative  method  using  a  region  of  interest  map  which  consists  of  either  only  vessels  or  specific 
brain  structures  (in  this  case  the  basal  ganglia)  plus  vessels,  a:  The  initial  xi-o(r)  map  without  noise  added  in  the  original  simulated 
images,  b:  xi=3 M  maP  °f  (a)  using  threshold  value  a  =  0.1.  c:  Similar  to  (b),  a  =  0.2.  d:  The  initial  x/=oM  map  with  noise  added  in  original 
images,  resulting  a  standard  deviation  of  0.025  radian  in  phase  images,  e:  x/=3(r)  map  of  (d)  using  a  =  0.1.  f:  The  associated  vessel 
map.  g:  The  x/=oM  map  in  the  coronal  plane  as  a  reference.  The  streaking  artifacts  are  clearly  shown  in  every  structure,  h:  The  \i= 3W 
maps  created  by  using  a  region  of  interest  map  which  consists  of  GR  SN,  RN,  PUT,  CN,  and  vessels,  i:  The  difference  image  of  (g)  and 
(h).  j:  The  initial  x/=o(r)  maP  in  the  transverse  plane  has  ‘false’  internal  capsule  (1C)  (pointed  by  an  arrow)  around  GP;  (k)  The  x/= 3(1)  map 
shows  no  “1C.”  This  matches  the  originally  simulated  model  (I).  No  noises  were  added  to  images  from  (g)  to  (I). 


from  a  uniform  region  inside  the  WM  decreased  from 
0.042  ppm  in  X;=oM  map  to  0.035  ppm  and  0.023  ppm 
in  the  Xr=3W  map  with  a  =  0.1  and  0.2,  respectively. 
The  baseline  susceptibilities  of  the  major  structures  are 
higher  with  SffARP  than  with  the  HP  filter.  The  iterative 
method  works  for  brain  structures  too  when  the  structure 
is  included  in  the  mask.  For  instance,  the  mean  suscepti¬ 
bility  values  of  the  GP  and  SN  have  been  increased  from 
0.155  ±  0.058  ppm  and  0.162  ±  0.067  ppm  in  the  Xr=oW 


map  to  0.163  ±  0.070  ppm  and  0.186  ±  0.083  ppm  in 
the  Xi— 3(f)  map,  from  the  dataset  processed  using 
SHARP.  The  result  after  HP  filtering  (Fig.  9e)  shows 
more  edge  artifacts  indicated  by  the  left  arrow  in  Fig.  9e. 
Much  of  this  error  was  reduced  by  the  iterative  method 
(Fig.  9f).  It  seems  that  the  iterative  method  compensated 
for  the  worse  first  guess  (Fig.  9e)  and  ended  up  with 
almost  the  same  result  (Fig.  9f,h)  as  having  started  with 
SHARP  (Fig.  9b, d)  from  the  image  perspective.  Since  a 
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FIG.  7.  The  plots  of  mean  susceptibility  values  inside  the  vessel 
(vein  of  Galen),  GP,  SN,  RN,  PUT,  and  CN  from  x/= 3W  maps.  The 
first  two  data  points  of  each  curve  is  the  value  inside  each  struc¬ 
ture  from  the  brain  model  and  the  X/=o(0  map,  respectively,  a: 
Xi= 3(r)  maps  generated  by  applying  a  region  of  interest  map  which 
consists  only  vessels  using  a  =  0.1,  0.15,  0.2,  0.25,  and  0.3, 
respectively,  b:  Xis(r)  maps  generated  by  applying  a  region  of 
interest  map  which  consists  of  the  GP,  SN,  RN,  PUT,  CN,  and 
vessels  using  a  =  0.1 . 


small  sized  HP  filter  cannot  remove  rapid  phase  wrap¬ 
ping  at  air-tissue  interfaces;  we  had  to  cut  out  the  region 
near  the  sinuses  in  the  phase  images. 


DISCUSSION 

In  this  article,  a  threshold-based  k-space/image  domain 
iterative  approach  has  been  presented.  Simulations  and 
in  vivo  results  show  that  the  ill-posed  problems  of 
streaking  artifacts  and  biases  in  the  estimates  of  suscepti¬ 
bilities  can  be  significantly  reduced.  The  replacement  of 
the  x(k)  values  near  the  singularities  by  XvmM-  which  is 
obtained  from  the  geometric  information  from  the  xM 
map  itself,  obviates  many  of  the  current  problems  seen 
in  the  TBSO  methods.  Since  XvmW  contains  little  streak¬ 
ing  artifacts  itself,  the  values  used  inside  the  thresholded 
regions  in  xffl  now  contain  no  artifact  either.  In  this 
sense,  we  obtain  an  almost  perfect  k-space  without  bad 
data  points  in  the  region  of  singularities.  This  explains 
why  this  method  converges  quickly  and  the  major 
improvement  is  in  the  first  iterative  step  (Fig.  5). 

The  proposed  iterative  approach  is  different  from  the 
other  threshold-based  methods  (11,15,18,19,24]  which 
fill  a  predefined  conical  region  using  a  constant,  zero  or 
1/thr  threshold  (11,15,24)  or  the  first-order  derivative  of 
g_1(k)  (19).  The  iterative  method  uses  full  geometry  in¬ 
formation  from  the  SM  (vessels  or  predefined  structures 
and  not  edge  information)  to  iteratively  change  k-space 
values  in  the  conical  region  using  the  forward  model. 
This  is  also  quite  different  than  other  currently  proposed 
solutions  (9,12,20,22).  Even  though  spatial  priors  such  as 
gradients  of  the  magnitude  are  used  (9,12,20,22),  in  those 
methods,  the  meaningful  values  of  the  singularity 
regions  in  k-space  are  obtained  through  solving  the  com¬ 
plex  cost  function  problem.  However,  the  iterative 
method  uses  priors  not  from  the  magnitude  image  but 
from  the  SM.  The  missing  data  in  the  singularity  regions 
are  obtained  through  iterating  back  and  forth  between 


0.02ppm 


-0.03ppm 


i 

FIG.  8.  Comparison  of  the  reconstructed  x/=3(r)  rnaps  using  (j)  accurate,  (k)  dilated,  and  (I)  eroded  vessel  maps.  Their  corresponding 
vessel  maps  and  the  enlarged  views  from  the  rectangular  regions  are  provided  in  (b)-(d)  and  (f)— (h).  (a)  and  (e)  The  initial  xi-o(r)  maps 
and  (i)  the  original  brain  model  as  references.  The  circle  in  the  midbrain  in  the  x(r)  maps  represents  the  RN  and  is  indicated  by  a  black 
arrow  in  (i).  Other  hyper-intense  regions  in  SMs  are  vessels. 
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FIG.  9.  Comparisons  of  SMs  using  SHARP  or  a  HP  filter  (FWHM  =  16  pixels)  to  remove  the  background  field.  The  iterative  method 
with  a  =  0.1  and  0.2  is  applied  after  the  background  is  removed,  (a)-(d)  and  (e)-(h)  are  results  after  the  application  of  SHARP  and  the 
HP  filter,  respectively,  (a)  and  (e)  the  initial  x/=oM  rnaps.  (b)  and  (f)  the  Xi=3(r)  maps  generated  from  the  iterative  method  with  a  =  0.1.  (c) 
and  (g)  the  differences  of  images  between  (a)  and  (b),  and  between  (e)  and  (f),  respectively.  These  two  images  show  the  successful 
reduction  of  the  streaking  artifacts,  (d)  and  (h)  the  x/= 3(1)  maps  generated  from  the  iterative  method  with  a  =  0.2.  The  range  of  the  gray¬ 
scale  bars  is  chosen  to  highlight  the  artifacts  in  the  images.  It  does  not  reflect  the  quantified  susceptibility  values  inside  veins. 


the  SMs  and  their  k-space.  The  advantage  of  cost  func¬ 
tion  approaches  is  that  they  do  not  need  to  predefine  the 
singularity  region  in  k-space  which  is  solved  by  the  opti¬ 
mization  process  automatically  (although  the  optimiza¬ 
tion  process  itself  is  usually  quite  time-consuming].  On 
the  other  hand,  the  iterative  method  is  the  most  time-ef¬ 
ficient.  It  is  fast  enough  to  reconstruct  SMs  for  a  512  x 
512  x  256  data  set  using  an  Intel  Core  i7  CPU  3.4  GHz 
processor  in  less  than  30  s,  since  in  practice  usually 
three  iterations  are  good  enough  to  generate  decent 
results. 

The  threshold  value  also  plays  a  key  role.  A  threshold 
value  of  0.1  is  a  reasonable  choice  since  a  lower  thresh¬ 
old  value  leads  to  an  increase  in  noise  and  a  higher 


threshold  value  leads  to  a  blurring  of  the  object  (Figs. 
6c,9d,h). 

It  is  known  that  the  ill-posedness  of  the  inverse  filter 
will  increase  the  noise  level  from  the  phase  to  the  SM. 
Based  on  both  simulations  and  real  data,  we  find  that 
there  is  a  factor  of  4  increase  in  noise  in  the  SM  relative 
to  the  original  phase  data.  This  result  and  the  fact  at  B0 
=  3  T,  Te  =  5  ms,  and  aXJ=0(r)  =  0.025  ppm  make  it  pos¬ 
sible  to  write  the  total  noise  in  the  background  region  in 
Xi=oM  as  0.025-4  (3/Bo)  (5/TE)  in  ppm.  The  noise  in 
Xj=3M  wiH  be  less  than  this  value  since  the  iterative 
method  will  reduce  streaking  artifact  in  SM. 

The  iterative  method  can  be  used  to  remove  streaking 
artifacts  associated  with  not  only  vessels  but  also  other 
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Table  1 

Ax  Measured  In  Vivo  in  Two  Veins  in  \j= o  if)  Maps  and  x/= 3  if) 
Maps  With  Different  Threshold  Values 


50= o  ( r ) 
map 

50=3  ( r ) 

map /a  =  0.1 

50=3  if) 
map /a  =  0.2 

VI  (SHARP) 

0.32  ±  0.07 

0.37  ±  0.08 

0.38  ±  0.09 

VI  (HP) 

0.24  ±  0.05 

0.28  ±  0.06 

0.28  ±  0.06 

V2  (SHARP) 

0.35  ±  0.04 

0.40  ±  0.05 

0.41  ±  0.05 

V2  (HP) 

0.25  ±  0.05 

0.31  ±  0.06 

0.30  ±  0.06 

Mean  and  standard  deviation  for  the  susceptibility  values  (in  ppm) 
of  two  veins  processed  using  SHARP  and  a  HP  filter  (FWHM  =  16 
pixels),  respectively,  were  chosen  from  the  0.5  mm  isotropic  reso¬ 
lution  data.  VI  and  V2  are  shown  in  Fig.  9.  The  susceptibility  val¬ 
ues  of  these  two  veins  have  been  increased  by  the  iterative 
method.  There  is  not  much  variation  of  the  susceptibility  value 
with  different  threshold  values. 


brain  structures  as  well.  Figure  6h  shows  a  reduction  in 
artifacts  associated  specifically  with  iron-rich  regions 
such  as  the  GP  and  CN. 

Accurately  extracting  vessels  from  the  Xr=o(r)  map  is 
critical  for  the  iterative  method  (Fig.  8).  In  this  study, 
vessels  were  segmented  directly  from  the  SM  (Fig.  2).  It 
may  also  possible  to  segment  veins  from  original  magni¬ 
tude  images  (9,12,20,22),  phase  images,  and/or  SWI 
images.  Extraction  of  accurate  anatomic  information  from 
phase  data  sometimes  is  difficult  since  phase  is  orienta¬ 
tion  dependent  and  phase  changes  are  generally  nonlo¬ 
cal.  SWI  images  work  better  for  an  anisotropic  dataset 
rather  than  an  isotropic  dataset  since  phase  cancellation 
is  needed  to  highlight  vessel  information.  Therefore,  we 
may  consider  combining  SMs  with  magnitude  images, 
phase  images,  and/or  SWI  images  together  to  segment 
the  veins,  since  different  types  of  images  can  compensate 
for  missing  information. 

The  iterative  method  appears  to  help  even  in  the  pres¬ 
ence  of  non-isotropic  resolution  with  partial  volume 
effects  and  to  a  minor  degree  when  an  HP  filter  is 
applied.  A  smaller  sized  HP  filter  would  be  better,  since 
a  larger  HP  filter  will  significantly  underestimate  the 
susceptibility  value  (Table  1).  SHARP  gave  us  better 
results  compared  with  the  HP  filter  (FWHM  =  16  pixels) 
(Fig.  9),  but  SHARP  requires  phase  unwrapping  which 
can  be  time  consuming  and  is  noise  dependent  (19). 
From  this  perspective,  an  HP  filter  has  the  advantage 
since  it  does  not  need  unwrapped  phase.  If  the  forward 
modeling  approach  of  Neelavalli  et  al.  (38)  can  be  used 
to  reduce  air/tissue  interface  fields,  then  it  may  be  possi¬ 
ble  to  use  a  small  size  HP  filter  (FWHM  =  8  pixels) 
which  may  provide  similar  results  to  SHARP. 

Severe  streaking  artifacts  associated  with  structures 
having  high  susceptibility  values  such  as  veins  can  lead 
to  major  changes  in  the  appearance  of  the  brain  struc¬ 
tures  with  low  susceptibility.  Practically,  the  susceptibil¬ 
ity  of  the  veins  is  a  factor  of  2.5-20  times  higher  than 
other  structures  in  the  brain.  Therefore,  even  a  10% 
streaking  artifact  can  overwhelm  the  information  in  the 
rest  of  the  brain  and  create  false  appearing  structures  as 
in  (Fig.  6j)  and  in  (Fig.  8e).  The  reduction  of  these  arti¬ 
facts  makes  a  dramatic  difference  in  the  ability  to  prop¬ 
erly  extract  the  susceptibility  of  other  tissues. 


In  conclusion,  both  simulations  and  human  studies 
have  demonstrated  that  the  proposed  iterative  approach 
can  dramatically  reduce  streaking  artifacts  and  improve 
the  accuracy  of  susceptibility  quantification  inside  the 
structures  of  interest  such  as  veins  or  other  brain  tissues. 
Given  its  relatively  fast  processing  time,  it  should  be 
possible  to  expand  its  use  into  more  daily  clinical  prac¬ 
tice.  With  the  improved  accuracy  of  the  susceptibility 
values  inside  veins,  this  method  could  be  used  poten¬ 
tially  to  improve  quantification  of  venous  oxygen  satura¬ 
tion  (18). 
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Purpose:  To  achieve  simultaneous  high-resolution  mag¬ 
netic  resonance  angiography  and  venography  (MRAV) 
imaging  in  terms  of  enhanced  time-of-flight  (TOF)  angiog¬ 
raphy  and  susceptibility-weighted  imaging  (SWI),  with  a 
clear  separation  of  arteries  and  veins. 

Materials  and  Methods:  A  new  flow  rephase /dephase 
interleaved  double-echo  sequence  was  introduced  that 
facilitates  multiple  types  of  imaging  contrast,  i.e.,  TOF. 
SWI,  and  dark  blood,  in  a  single  acquisition.  A  nonlinear 
subtraction  (NLS)  method  is  proposed  and  assessed  to 
maximally  enhance  angiography  contrast  with  minimum 
venous  contamination. 

Results:  Fully  flow  rephased  TOF  MRA  and  SWI  MRV 
data  were  acquired  simultaneously,  along  with  an  extra 
flow  dephased  dark  blood  image  for  angiography 
enhancement  calculation.  Compared  to  linear  subtraction 
methods,  the  proposed  NLS  method  was  able  to  enhance 
angiography  contrast  while  effectively  eliminating  vein- 
tissue  contrast.  The  NLS  method  even  outperformed 
low-dose  contrast-enhanced  MRA  for  a  clean,  enhanced 
angiography  map. 

Conclusion:  Using  the  proposed  interleaved  double-echo 
sequence  along  with  the  NLS  postprocessing  method,  one 
can  simultaneously  obtain  both  high-quality  SWI  and  sig¬ 
nificantly  enhanced  TOF  MRA  with  clear  separation  of  ar¬ 
terial  and  venous  maps. 
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UNDERSTANDING  THE  CEREBRAL  vascular  system 
morphologically  and  functionally  is  very  important  for 
both  neuroscientific  studies  and  clinical  applications. 
Numerous  clinical  problems  such  as  arteriovenous  mal¬ 
formation  (AVM)  (1),  tumor  evaluation  (2-3),  and  trau¬ 
matic  brain  injury  (TBI)  (4,5)  require  detailed  vascular 
information  for  best  diagnostic  interpretation.  Due  to 
the  vast  differences  in  hemodynamic  characteristics, 
pathological  behavior,  and  the  underlying  magnetic 
resonance  imaging  (MRI)  principles  of  vessel-tissue  con¬ 
trast,  arteries  and  veins  are  usually  imaged  with  sepa¬ 
rate  protocols,  leading  to  increased  scanning  time  and 
motion-related  misregistrations.  Recently,  there  has 
been  a  focus  on  single-scan  magnetic  resonance  angi¬ 
ography  and  venography  (MRAV)  imaging  (6-9)  in  order 
to  obtain  a  more  complete  perspective  of  the  vascula¬ 
ture.  Collecting  MRA  and  MRV  in  a  single  acquisition 
gets  rid  of  any  potential  registration  problems  and 
enables  a  precise  review  and  a  complete  picture  of  the 
relationship  between  the  two  vasculature  networks. 

However,  further  applications  of  these  MRAV  meth¬ 
ods  are  hindered  by  two  major,  yet  conflicting,  objec¬ 
tives:  enhancing  vessel-tissue  contrast  and  distin¬ 
guishing  arteries  and  veins.  In  MRA,  bright-blood  (BB) 
strategies  such  as  time-of-flight  (TOF)  and  contrast- 
enhanced  MRA  (CE  MRA)  are  used  to  obtain  high 
blood  signal  while  reducing  tissue  signal  as  much  as 
possible.  MRV,  on  the  other  hand,  can  be  most  reli¬ 
ably  done  with  dark-blood  (DB)  methods  such  as  sus¬ 
ceptibility-weighted  imaging  (SWI)  (10-12)  to  selec¬ 
tively  suppress  venous  signal  while  maintaining  high 
tissue  signal  for  better  vein-tissue  contrast  or  arterial 
saturation.  Since  the  tissue  signal  is  unwanted  in 
MRA  scan  but  is  required  for  optimal  venous  contrast 
in  SWI,  obtaining  both  simultaneously  is  a  challenge, 
regardless  of  whether  using  single-echo  (7),  double¬ 
echo  (6),  or  multiecho  (9)  strategies.  Gadolinium  (Gd)- 
based  contrast  agents  can  be  used  as  either  a  BB  or 
DB  (13)  method  that  can  enhance  both  arteries  and 
veins  at  the  same  time  due  to  the  global  vascular  T !  or 
T2  shortening  effects;  thus,  its  use  will  make  it  more 
difficult  to  separate  arteries  and  veins  in  CE  MRA. 

Currently  the  gold  standard  of  angiography  is  intra¬ 
arterial  digital  subtraction  angiography  (IADSA)  (14,15). 
The  basic  principle  is  to  create  two  different  levels  of 
blood  signal  while  keeping  tissue  signal  the  same  by 
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Figure  1.  Diagram  of  the  FR/FD  inter¬ 
leaved  double-echo  sequence.  Echon, 
Echo12,  and  Echo2i  are  flow  rephased 
(FR)  using  first-order  moment  nulling 
gradients  dynamically  calculated  with 
the  phase  encoding  gradients  and  TEs 
taken  into  account,  and  Echo12  is  flow 
dephased  (FD)  by  bipolar  gradients 
with  a  low  VENC  value. 


administrating  an  iodinated  contrast  agent  into  the 
intravascular  space,  and  generate  a  pure  vessel  map 
simply  using  subtraction.  Similar  attempts  have  been 
made  to  enhance  MRA  contrast  by  using  a  rephase/ 
dephase  strategy  to  create  high/low  blood  signals  while 
maintaining  the  tissue  signal  (16-20).  In  the  recent 
work  by  Kimura  et  al  (16),  a  double-echo  sequence  col¬ 
lecting  flow  rephased/dephased  images,  called  hybrid 
of  opposite  contrast  MRA  (HOP-MRA),  was  developed 
and  a  linearly  weighted  subtraction  was  proposed  to 
increase  the  vessel-background  contrast  by  minimizing 
the  subtracted  tissue  signals.  Although  capable  of 
enhancing  arteries,  the  major  issue  of  linear  subtrac¬ 
tion,  weighted  or  not,  is  that  it  will  introduce  heavy  ve¬ 
nous  contamination  such  that  veins  also  appear  bright 
in  the  subtracted  images.  As  a  result,  it  is  difficult  to 
distinguish  between  arteries  and  veins,  especially  for 
smaller  vessels.  The  venous  contamination  can  be  pos¬ 
sibly  reduced  if  both  flow  rephased  and  dephased 
images  are  acquired  with  long  (preferably  equal)  TEs,  so 
that  veins  are  already  suppressed  because  of  Tj  decay. 
With  equal  TE,  the  rephase/ dephase  subtraction  can 
reduce  the  tissue  signal-to-noise  level,  and  can  better 
reveal  the  arteries  after  processing  such  as  maximum 
intensity  projection  (MIP)  or  vessel  tracking.  However, 
flow  compensation  at  long  TE  will  be  less  efficient,  espe¬ 
cially  for  fast  flows,  leading  to  artifacts  and  signal  non¬ 
uniformity  in  large  arteries  such  as  carotid  arteries  (21). 

To  tackle  this  challenge,  we  propose  using  a  flow 
rephase /dephase  interleaved  double-echo  gradient 
recalled  echo  (GRE)  sequence  capable  of  MRAV  imag¬ 
ing  and  at  the  same  time  acquiring  both  flow 
rephrased  and  dephased  images.  A  nonlinear  subtrac¬ 
tion  processing  method  is  then  proposed  for  selective 
MRA  enhancement  utilizing  the  flow  rephrased  and 
dephased  images.  The  performance  of  both  nonlinear 
and  linear  subtraction  methods  will  be  analyzed  and 
compared  for  the  arteiy-tissue  contrast  theoretically 
and  experimentally. 

MATERIALS  AND  METHODS 
Sequence  Design 

In  order  to  acquire  MRAV  data  as  well  as  flow 
dephased  data  for  MRA  enhancement  (16),  we  devel¬ 
oped  a  flow  rephase/dephase  interleaved  dual-echo 
GRE  sequence  as  shown  in  Fig.  1.  In  the  first  TR, 


both  echoes  were  fully  flow  rephased  (FR,  or  flow  com¬ 
pensated)  on  all  three  axes,  while  in  the  second  TR 
the  second  echo  was  flow  dephased  (FD)  by  a  series  of 
bipolar  gradient  pairs. 

Unlike  Du  and  Jin’s  (6)  method  in  which  flow 
rephasing  for  the  second  echo  was  achieved  with  a 
fly-back  gradient  in  the  readout  direction,  we  imple¬ 
mented  an  active  scheme  for  flow  rephasing  to  the 
second  echo  on  all  three  axes  by  dynamically  calculat¬ 
ing  the  gradients  for  first-order  moment  nulling  on 
each  axis,  taking  into  account  TE  as  well  as  phase 
encoding  or  readout  gradients  (22).  In  this  way,  one 
can  independently  set  the  two  echoes  for  their  TE, 
bandwidth,  and  asymmetry  without  affecting  the  flow 
rephasing  performance.  Flow  dephasing  for  the  sec¬ 
ond  echo  was  done  by  using  bipolar  gradients  of  very 
low  velocity  encoding  (VENC)  value. 

With  this  sequence,  four  images  will  be  generated 
with  a  single  scan:  two  FR  images  with  a  short  TEi  for 
normal  TOF-MRA  (Echon  and  Echo2i),  one  FR  image 
with  a  long  TE2  for  SWI  (Echo12),  and  one  FD  image 
(Echo22)  also  with  TE2.  Since  excellent  SWI  results 
can  be  reliably  generated  from  the  Echo12  data  (6,10), 
the  following  analysis  will  focus  on  enhancing  MRA 
contrast  and  separating  venous  and  arterial  signal  by 
assessing  and  comparing  nonlinear  subtraction  and 
linear  subtraction  methods. 

Linear  Subtraction  (LS)  vs.  Nonlinear 
Subtraction  (NLS) 

Since  the  contrast  characteristics  of  linear  subtraction 
have  been  extensively  discussed  in  previous  studies, 
interested  readers  are  referred  to  Kimura  et  al' s  article 
for  details  (16).  Here,  we  will  further  assess  the  effect 
of  different  weighting  factors  in  the  subtraction  and 
their  effects  on  vessel-tissue  contrast. 

Denoting  the  signal  intensity  under  two  different 
imaging  conditions  (eg,  FR  and  FD)  as  S  and  S’,  and 
assuming  the  same  noise  standard  deviation  (SD)  a 
for  both  signal  levels,  one  can  calculate  the  signal-to- 
noise  ratio  (SNR)  of  the  LS  image  as: 

SNRis  =  A5b  =  =  [1] 

oV  1  +  a2  oV  1  +  a2  a\/l  +  a2 

where  a  is  the  weighting  factor  for  the  subtraction, 
and  |3  =  S’/S  is  defined  as  the  ratio  between  the  two 
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Table  1 

Estimated  Signal  Ranges  of  Arteries,  Veins,  and  Tissue  in  FR  and 
FD  Images 


Artery 

Vein 

Tissue 

S  (FR) 

High 

Medium 

Medium 

S’  (FD) 

Low 

Low 

Medium  but  S’  <  S 

signal  levels.  The  CNR  between  artery  and  tissue,  i.e., 
CNRat  is,  can  be  written  as: 


CNRatls  =  ^ 


3.)  S') 

<J\/ 1  +  a2 


[2] 


The  subscripts  “a”  and  “t”  denote  artery  and  tissue, 
respectively.  The  vein-tissue  CNRvt.is  can  be  written  in 
the  same  form.  On  the  other  hand,  if  the  signal  is 
nonlinearly  weighted  in  a  self-weighted  manner  (i.e., 
squared  on  a  voxelwise  basis)  prior  to  subtraction  so 
that  ASnis  =  S2  —  aS’2,  the  noise  after  NLS  can  be 
approximated  by: 


Unis  =  2  \J  (S2  +  a2S'2){7 


and  SNRnls  can  be  written  as: 


SNRnis  = 


S2  -  aS'2 


Onls 


S(l-aP2) 

2  u\J  1  +  a2p2 


[3] 

[4] 


Comparing  with  Eq.  [1],  the  (3  factor  now  appearing 
in  the  denominator  in  Eq.  [4]  is  the  major  difference 
between  LS  and  NLS.  The  value  of  (3  depends  on  the 
type  of  tissue,  i.e.,  whether  they  are  arteries,  veins,  or 
parenchymal  tissues,  and  it  varies  according  to  their 
respective  rephase/dephase  signals  and  partial 
volume  effects. 

The  expression  for  CNRnis  between  vessels  and 
tissues  is  much  more  complicated  than  the  LS  version 
shown  in  Eq.  [2].  From  a  qualitative  point  of  view,  S’a 
and  S’v  will  be  very  low  or  even  nulled  as  a  combined 
result  of  low  VENC  gradients  along  with  high  flow  ve¬ 
locity  (for  arteries)  or  T *2  decay  (for  veins),  while  Sa  will 
be  higher  than  Sv  thanks  to  greater  TOF  inflow.  For 
static  parenchymal  tissues,  on  the  other  hand,  the 
dephasing  gradients  will  not  significantly  change  the 
signal  as  its  equivalent  diffusion  weighting  effects  are 
negligible  (see  next  session).  However,  since  in  our 
sequence  S'  is  collected  with  longer  TE  so  that  it  will 
decay  (relative  to  S  collected  at  short  TE)  by  ~20% 
with  a  typical  ATE  of  15  msec  and  a  tissue  T2  of  80 
msec  at  3T,  so  that  (3  <  1  remains  true  for  tissues. 
The  summary  of  the  above  analysis  is  listed  in  Table 
1 .  Generally  speaking,  (3  will  be  of  low  value  for  both 
arteries  and  veins  (especially  for  arteries  due  to  the 
high  Sa),  but  is  expected  to  be  around  0.8  for  tissues 
mainly  due  to  Tj  decay. 

To  confirm  this,  we  performed  Monte  Carlo  simula¬ 
tions  according  to  Eqs.  [1-4]  to  evaluate  the  signal 
change  ASis  and  ASnis,  noise  SD,  and  SNR  for  all  S> 
S’  (i.e.,  (3  <  1)  combinations.  S  and  S'  ranged  from  1 
to  300  and  were  added  with  Gaussian  noise  (ct  =  10), 
yielding  a  maximum  SNR  of  about  30  in  the  original 


signals.  Also,  simulation  on  the  effect  of  weighting  fac¬ 
tor  a  to  ASis  and  ASn]s  are  also  performed  to  estimate 
and  predict  its  optimal  value. 


Data  Acquisition  and  Processing 

Six  healthy  volunteers  (two  males  and  four  females, 
age  27  ±  4  years)  were  recruited  for  the  study  and 
written  consents  approved  by  the  local  Institutional 
Review  Board  were  obtained  from  each  volunteer  prior 
to  the  scans.  All  scans  were  performed  on  a  Siemens 
3T  Verio  system  (Siemens  Healthcare,  Erlangen,  Ger¬ 
many)  with  a  product  32-channel  head  coil.  The  scan¬ 
ning  parameters  of  the  FR/FD  interleaved  double¬ 
echo  GRE  sequence  were:  TEi/TE2/TR  =  6.2/22/28 
msec,  flip  angle  (FA)  =  15°,  bandwidth  =  260  Hz/ 
pixel,  voxel  size  =  0.5  x  0.5  x  1.0  mm3,  slice  number 
=  48,  and  total  scan  time  was  6  minutes  40  seconds. 
The  bipolar  gradients  for  flow  dephasing  were  4  msec 
in  duration  and  24  mT/m  in  amplitude,  leading  to  a 
VENC  value  of  about  1.46  cm/s.  As  a  side  note,  previ¬ 
ous  literature  on  similar  studies  used  the  b  value  (16) 
rather  than  the  VENC  value  for  describing  the 
dephasing  effects.  However,  the  b  values  are  actually 
too  small  to  introduce  significant  difference  observed 
in  the  vessels  (eg,  b  =  2s/mm2  reported  previously 
(16)  and  b  =  1.8  s/mm2  in  our  case),  thus  the  blood 
signal  loss  in  the  FD  image  is  primarily  due  to  flow 
dephasing  effects.  One  of  the  subjects  also  addition¬ 
ally  underwent  low-dose  contrast-enhanced  MRA 
scans  with  an  injection  of  5  cc  Magnevist  (gadopente- 
tate  dimeglumine) ,  and  single  echo  flow  rephased 
GRE  sequence  was  scanned  using  a  matching  proto¬ 
col  (with  TE/TR  =  6.2/28  msec  and  FA  =  15°)  to  col¬ 
lect  CE  TOF  MRA  images  as  comparison  to  the  sub¬ 
traction  results.  GRAPPA  with  an  acceleration  factor 
of  2  was  used  for  all  scans. 

SWI  results  were  calculated  from  the  flow  rephased 
Echo12  data,  using  a  highpass  filter  of  64  x  64  on  the 
phase  image  and  a  four -time  phase  mask  multiplica¬ 
tion  (i.e.,  m  =  4)  (10).  mIP  images  over  a  12-mm  slab 
were  created  to  better  visualize  the  veins  as  continu¬ 
ous  structures. 

For  both  LS  and  NLS  calculations,  the  two  identical 
short  TE  FR  images  (i.e.,  Echon  and  Echo21)  were 
averaged  to  improve  SNR  especially  for  arteries,  and 
the  averaged  image  was  used  as  S.  Taking  the  FD 
image  as  S’,  LS  and  NLS  results  were  then  calculated 
with  a  =  0.75,  1.00,  and  1.25,  respectively.  MIP 
images  over  a  12-mm  slab  were  also  created  for  all 
final  results  to  display  the  arteries.  The  skull  signal 
was  removed  for  better  visualization  using  SPIN  (Sig¬ 
nal  Processing  In  NMR,  Detroit,  Michigan),  and  a 
brain-only  binary  mask  was  obtained  and  applied  to 
all  processed  data  described  above  (but  prior  to  inten¬ 
sity  projections)  to  remove  the  skull. 

The  artery-tissue  CNR  of  different  artery  sizes  was 
also  determined  to  evaluate  the  artery  enhancing  per¬ 
formance  on  various  vessel  sizes.  The  contrast  was 
determined  by  subtracting  the  mean  signal  intensity  of 
the  surrounding  tissue  from  the  signal  maxima  of  the 
artery  segment,  and  noise  was  taken  as  the  SD  of  the 
tissue  signal.  The  artery  diameter  was  estimated  as  the 
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Figure  2.  Simulation  results  of  the  subtraction  signal  (ASis  in  (a)  and  ASnis  in  (d)),  noise  SD  (b,e),  and  the  resulting  SNR  (c,f)  of 
LS  (a-c)  and  NLS  (d-f).  The  simulation  was  performed  for  all  S  >  S’  (i.e.,  (3  <  1)  combinations,  and  a  =  1  and  cr  =  10  was  used 
for  calculations  here.  Higher  temperature  color  corresponds  to  higher  subtraction  signal  (a,d),  or  higher  noise  SD  (b,e),  or 
higher  SNR  of  the  subtracted  image  (c,f).  The  three  ROIs  indicate  the  estimated  regions  for  arteries  (“A"),  veins  (“V”)  and  tissues 
(*T”)  signals.  The  color  scales  of  NLS  were  scaled  to  match  the  range  of  LS  for  the  columns  of  subtracted  signal  and  noise  SD. 


half-wldth-maximum-height  of  its  profile  and  was  sorted 
into  several  groups:  less  than  1  mm,  1-1.5  mm,  1.5-2 
mm,  and  2-2.5  mm.  For  each  group,  live  vessels  were 
selected  and  the  mean  and  SD  of  CNR  were  evaluated. 


RESULTS 

The  simulation  comparing  LS  and  NLS  are  shown  in 
Fig.  2,  where  the  estimated  signal  regions  for  arteries, 
veins,  and  tissues  are  indicated  according  to  Table  1. 
The  simulation  confirmed  that  the  introduction  of 
factor  |3  into  the  denominator  in  Eq.  [4]  modifies  the 
nonlinear  subtraction  outcomes  between  rephased  and 
dephased  signals.  Specifically,  the  gradient  of  ASnis  is 
shifted  towards  high  S  when  S’  values  are  low,  but  the 
shift  is  much  less  significant  for  medium  to  high  S’  val¬ 
ues  (Fig.  2d).  For  LS,  ASis  of  the  veins  is  higher  than 
that  of  the  tissues,  thus  showing  a  positive  vein- tissue 
contrast  (Fig.  2a).  However,  venous  ASnis  is  located  at 
the  same  gradient  level  with  the  tissues,  effectively 
reducing  or  even  eliminating  the  vein-tissue  contrast. 
Meanwhile,  both  ASis  and  ASnis  of  the  arteries  remain 
very  high  (hot  colors),  maintaining  a  high  contrast  to 
the  tissues  (Fig.  2d).  In  light  of  the  above  results,  it  is 
predicted  that  both  veins  and  arteries  will  be  shown  in 
the  LS  approach,  while  predominantly  only  arteries 
will  be  visible  in  the  NLS  approach. 

The  weighting  factor  a  was  introduced  in  order  to 
best  null  the  subtracted  tissue  signal  (16).  The  effect 


of  a  is  to  “tilt”  the  gradient  of  ASis  or  ASn|s,  as  demon¬ 
strated  in  Fig.  3,  where  simulations  using  a  =  0.5, 
1.0,  and  1.5  were  compared.  When  a  higher  a  value  is 
used,  ASis  will  reduce  more  at  high  S’  values.  For 
practical  considerations,  since  S’  of  tissues  is  higher 
than  both  arteries  and  veins  in  the  dephased  images, 
one  can  expect  to  see  a  more  decreased  tissue  signal 
in  the  LS  images  with  higher  a  value.  This  is  consist¬ 
ent  with  the  suggested  optimal  value  of  1.5  for  a  (16). 
However,  the  difference  in  ASis  between  veins  and 
arteries  will  also  become  greater  with  high  a  value 
due  to  the  reduced  tissue  signal  (Fig.  3a-c),  suggest¬ 
ing  an  even  increased  vein-tissue  contrast  in  the  LS 
results.  For  the  nonlinear  approach,  ASnis  of  both 
veins  and  tissues  only  vary  slightly  between  a  values, 
and  the  vein-tissue  contrast  in  the  NLS  results  is  min¬ 
imal  and  insensitive  to  the  weighting  factor  (Fig.  3e,f). 

The  single  slice  images  of  the  unprocessed  FR  and 
FD  data  along  with  the  LS/NLS  processed  data  of  one 
representative  subject  are  shown  in  Fig.  4,  while  the 
corresponding  MIP  or  mIP  results  are  shown  in  Fig.  5. 
These  images  show  several  representative  arteries  and 
veins  that  can  be  used  as  a  demonstration  of  our 
approach.  The  profiles  of  a  pair  of  neighboring  artery 
and  vein,  normalized  by  the  respective  maximum  sig¬ 
nal,  were  extracted  as  indicated  by  the  thin  black  line 
in  Fig.  4a  and  shown  in  Fig.  6,  and  the  vessel-tissue 
contrast  and  CNR  value  of  these  vessels  are  shown 
in  Table  2.  The  contrast  was  calculated  as  the  signal 
difference  between  vessel  center  and  surrounding 


Noncontrast  MRAV  Imaging 


5 


U) 


Figure  3.  Simulation  results  of  the  subtraction  signal  ASis  (a-c)  and  ASnia  (d-f),  with  a  =  0.5  (a,d),  1  (b,e),  and  1.5  (c,f), 
respectively.  The  color  gradients  of  NLS  were  scaled  to  match  the  range  of  LS. 


tissues  divided  by  the  tissue  signal  mean,  and  the 
CNR  was  calculated  as  the  signal  difference  between 
vessel  center  and  surrounding  tissues  divided  by  the 
standard  deviation  of  tissue  signal.  The  (3  value  was 
about  0.07  and  0.27  for  the  arteiy  and  the  vein  (cal¬ 
culated  at  the  vessel  profile  center),  respectively,  and 
was  about  0.79  for  the  surrounding  tissues,  which  is 
in  agreement  with  the  above  theoretical  analysis.  The 
comparison  between  LS,  NLS,  and  low-dose  CE  TOF 
MRA  are  similarly  shown  in  Fig.  7,  and  the  artery-tis¬ 
sue  CNR  of  different  arteiy  sizes  are  shown  in  Fig.  8. 


DISCUSSION 

In  this  study  we  proposed  an  MRAV  imaging  method 
with  selective  MRA  enhancement.  Specifically,  we 
designed  an  FR/FD  interleaved  double-echo  GRE 
sequence  which,  for  the  first  time,  offers  the  capacity 
of  obtaining  both  fully  flow  rephased  SWI  data  for 
venography  and  TOF  MRA  with  selectively  enhanced 
angiography  and  minimal  venous  contamination.  Spe¬ 
cifically,  a  dynamic  full  flow  compensation  scheme 
and  a  nonlinear  subtraction  method  were  proposed  to 
achieve  the  above  objectives. 

For  noncontrast-enhanced  MRA,  the  3D  TOF 
method  is  most  widely  used  for  clinical  exams  for  its 
potential  to  image  small,  tortuous  arteries  noninva- 
sively.  TOF  MRA  creates  a  positive  arteiy-tissue  con¬ 
trast  as  a  result  of  the  in-flowing  arterial  blood  being 
much  less  saturated  than  the  static  tissues.  However, 
the  in-flowing  blood  will  gradually  become  saturated 
as  it  passes  through  the  imaging  slab,  losing  the 


contrast  especially  for  smaller  vessels  with  slow  flow. 
There  are  several  means  to  increase  the  contrast  for 
small  arteries,  including  using  magnetization  transfer 
contrast  (MTC)  radiofrequency  (RF)  pulse  to  suppress 
tissue  signal  (23,24),  multiple  thinner  slabs  to  reduce 
saturation  effect  (25),  TONE  pulse  for  spatially  varying 
flip  angle  excitation  (26),  or  acquiring  an  additional 
flow  dephased  image  and  then  subtract  it  from  the 
flow  rephased  image  (16-20).  The  subtraction  method 
is  of  particular  interest  as  it  can  greatly  reduce  or 
even  null  the  background  tissue  signal  while  retaining 
high  vessel  signal.  Furthermore,  a  common  dilemma 
for  existing  MRAV  methods  is  the  conflicting  require¬ 
ment  on  tissue  signal,  that  TOF  MRA  requires  low  tis¬ 
sue  signal  while  SWI-based  MRV  requires  high  tissue 
signal.  This  now  can  be  perfectly  addressed  with  our 
interleaved  double-echo  MRAV  sequence,  as  one  can 
have  sufficient  tissue  signal  for  SWI  using  a  proper 
flip  angle  and  have  reduced  tissue  signal  for  MRA  af¬ 
ter  subtraction. 

The  major  drawback  of  the  LS  method  is  that  as  the 
venous  blood  is  also  suppressed  by  the  dephasing  gra¬ 
dients,  vein-tissue  contrast  will  become  positive  and 
confounds  the  interpretation  of  the  MRA  results.  Previ¬ 
ous  studies  mainly  focused  on  maximizing  the  artery- 
tissue  contrast,  but  failed  to  address  such  venous  con¬ 
tamination  issue  (16,17).  Compared  to  LS  methods, 
the  proposed  NLS  method  is  proven  to  be  able  to 
enhance  the  artery-tissue  contrast  while  suppressing 
the  venous-tissue  contrast  (Figs.  4-7).  Self- weighted 
NLS  reduces  the  vein-tissue  contrast  by  additionally 
taking  into  account  the  signal  intensity  of  both 
rephased  and  dephased  states,  i.e.,  tissues  have  low 
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Figure  4.  Top  row:  single  slice  images 
of  (a)  TOF  images  from  the  averaged  FR 
echoes,  (b)  SWI  from  Echoi2,  (c)  FD 
dark  blood  data  from  Echo22.  Middle 
row:  LS  results  using  (d)  a  =  0.75,  (e)  a 
=  1.00,  and  (f)  a  =  1.25.  Bottom  row: 
NLS  results  using  (g)  a  =  0.75,  (h)  a  = 
1.00,  and  (i)  a  =  1.25.  Black  and  white 
arrows  indicate  two  veins  and  one  ar¬ 
tery  of  interest;  the  thin  black  line  indi¬ 
cates  the  location  at  which  the  vessel 
profiles  are  drawn  for  Fig.  6.  The  dis¬ 
play  window  was  the  same  among  each 
row,  except  for  the  SWI  image,  which 
was  individually  adjusted. 


S  —  S'  but  high  S  +  S’  while  veins  have  high  S-S’  but 
low  S+S’,  and  therefore  both  have  similar  S2  —  S’2. 

Apart  from  the  linearity  or  nonlinearity,  the  weight¬ 
ing  factor  a  in  the  subtraction  also  alters  the  vessel- 
tissue  contrast  by  changing  the  AS  distribution,  as 
demonstrated  in  Fig.  3.  With  a  greater  a  factor,  the  LS 
method  indeed  was  able  to  enhance  arteiy-tissue  con¬ 
trast  by  reducing  tissue  signal,  but  also  further 
enhanced  the  vein-tissue  contrast  (Fig.  6b).  These 
enhanced  veins  may  be  misinterpreted  as  small 
arteries  (Fig.  5)  that  have  medium  to  low  contrast  due 
to  partial  volume  effects  or  steady-state  saturation 
effects.  The  vein-tissue  contrast  is  minimized  in  NLS 
due  to  the  fact  that  both  veins  and  tissues  are  of  simi¬ 
lar  ASnis  level,  and  the  slight  variation  of  tissue  signal 
between  different  a  values  does  not  significantly 
change  their  contrast.  Although  using  a  =  0.75  in  LS 
can  also  lead  to  minimal  vein-tissue  contrast,  the 
artery-tissue  CNR  is  reduced  by  over  23%  compared 
to  a  =  1.25  (Fig.  6b  and  Table  2).  Note  that  there  was 
a  large  increase  in  CNR  of  LS  with  a  =  1.25.  This  was 
because  the  tissue  signal  was  completely  subtracted 


out  since  «p~l ,  thus  the  noise  distribution  changed 
from  Gaussian  to  Raleigh  and  resulted  in  a  smaller 
standard  deviation.  For  NLS,  in  contrast,  such  varia¬ 
tion  in  CNR  is  less  than  3%  between  all  a  values,  and 
a  around  1  is  seen  to  sufficiently  eliminate  the  vein- 
tissue  contrast  (Fig.  6c  and  Table  2).  We  observed 
such  insensitivity  of  a  in  NLS  results  in  all  subjects. 
This  means  that  one  simply  needs  to  square  the 
rephased/dephased  images  and  subtract  them  in 
order  to  obtain  clean,  enhanced  angiography. 

To  better  demonstrate  the  MRA  enhancing  effects  of 
LS  and  NLS,  we  also  collected  low-dose  CE  TOF  MRA 
data  for  comparison.  By  using  low-dose  contrast 
agent  (5  cc),  low  flip  angle  (15°),  and  long  TR  (28 
msec),  the  Tj  shortening  effect  of  the  contrast  agent  is 
mainly  utilized  to  minimize  saturation  of  the  blood 
signal  so  as  to  maximize  TOF  effects,  rather  than 
using  the  steady-state  Tj  enhancement,  as  in  full- 
dose  CE  MRA.  Although  Ti  weightings  may  still  con¬ 
tribute  more  or  less  to  the  enhanced  blood  signal, 
such  an  effect  is  not  strong,  as  evidenced  by  the  lack 
of  most  slow-flowing  veins  in  Fig.  7d  and  the  identical 
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Figure  5.  Corresponding  MIP  (or  mIP) 
images  of  Fig.  4.  The  white  and  black 
arrows  in  (d)  indicate  an  artery  and  a 
vein,  respectively.  All  results  shown 
were  projected  over  a  1-2-mm  slab. 


artery-tissue  contrast  between  normal  TOF  and  5  cc 
data  in  Fig.  7e.  Were  the  'I',  shortening  effect  strong, 
veins  would  also  have  become  visible  and  artery-tis¬ 
sue  contrast  would  be  higher.  By  this  means,  one  can 
establish  an  upper  baseline  for  maximal  TOF  effects 
in  arteries  to  compare  the  performance  of  LS  and  NLS 


with.  As  shown  in  Fig.  7,  arteries  are  enhanced  in  all 
three  results  being  compared.  Most  noticeably,  NLS 
seems  to  even  outperform  the  CE  TOF  data  by  its 
higher  artery-tissue  contrast  (Figs.  7e,  8),  extra  reve¬ 
lation  of  some  small  arteries  (e.g.,  near  the  frontal 
part),  uniform  background  signal,  and  the  lack  of 


Figure  6.  Normalized  vessel  profiles  of  a  vein  and  an  artery  from  (a)  unprocessed  FR  (TOF)  and  FD  (dark  blood)  data,  (b)  LS 
processed  data,  (c)  NLS  processed  data,  and  (d)  comparison  between  LS  and  NLS.  The  arrows  indicate  the  center  of  the  vein, 
which  will  not  be  seen  in  some  cases  due  to  reduced  contrast  with  surrounding  tissues.  Quantification  of  contrast  and  CNR 
are  shown  in  Table  2.  The  diameter  of  this  artery  was  ~2  mm. 
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Table  2 

Vessel-Tissue  Contrast  of  the  Artery  and  the  Vein  Whose  Profile  Was  Indicated  by  the  Black  Line  in  Fig. 4a 


Artery 

Vein 

Contrast  (%) 

CNR 

Contrast  (%) 

CNR 

TOF(FR) 

LS  (or=0.75/1 .00/1 .25) 

NLS  ( a=0.75/1 .00/1 .25) 

28.3 

206.1/453.8/2493.3 

255.4/389.4/680.8 

5.5 

17.8/17.8/23.2 

15.5/15.6/15.2 

-37.2 

-3.5/37.4/366.4 

-35.2/-20.0/12.6 

-7.2 

-0.3/1 .5/3.4 
-2.1/-0.8/0.3 

veins  or  the  contrast-enhanced  choroid  plexus.  LS 
offers  similar  artery-tissue  contrast  with  NLS  but  also 
lots  of  venous  signal,  making  it  very  much  like  the 
routine  full-dose  CE  MRA  data  (data  not  shown).  The 
comparison  confirms  that  the  proposed  self-weighted 
NLS  method  is  potentially  capable  of  revealing  the 


fel4  & 


(c)  (d) 


Figure  7.  MIP  images  of  (a)  TOF  MRA.  (b)  LS,  (c)  NLS,  and 
(d)  low-dose  CE  MRA  with  5  cc  gadopentetate  dimeglumine. 
(a-c)  Collected  with  a  single  scan  of  the  FR/FD  interleaved 
double-echo  sequence  before  the  use  of  contrast  agent,  and 
a  was  1  for  both  LS  and  NLS  calculation,  (e)  Comparison  of 
the  normalized  profiles  of  one  representative  artery  as  shown 
in  (a),  extracted  from  the  corresponding  single  slice  images. 


maximal  arterial  TOF  effect  without  the  use  of  con¬ 
trast  agent.  One  thing  worth  mentioning  is  that  due  to 
its  fast  flow,  the  superior  sagittal  sinus  would  also 
appear  bright  in  all  methods  compared  (Fig.  7),  albeit 
it  can  be  easily  distinguished.  Inferior  sagittal  sinus, 
however,  is  only  seen  with  contrast  agent  (Fig.  7d)  but 
not  in  TOF,  LS,  or  NLS  results. 

Due  to  partial  volume  effects  and  blood  flow  veloc¬ 
ity,  it  is  expected  that  the  CNRat  may  vary  according 
to  artery  sizes.  Moreover,  because  squaring  the  image 
in  NLS  will  amplify  the  noise  by  twice  the  original  sig¬ 
nal  intensity,  SNRnls  will  have  lower  value  than  SNRls 
(Fig.  2f  vs.  2c)  on  a  voxelwise  basis.  According  to  Fig. 
8,  CNRat  becomes  greater  for  larger  vessels  with  all 
methods,  which  is  expected.  Paired  t-test  results  indi¬ 
cate  that  CNRat  of  NLS  is  significantly  greater  than  LS 
for  arteries  larger  than  1.5  mm  (two-tailed  P  =  0.046 
and  0.035  for  1. 5-2.0  mm  and  2. 0-2. 5  mm  groups, 
respectively),  while  both  have  similar  CNRat  for  those 
smaller  than  1.5  mm  (P  =  0.126  and  0.395  for  <lmm 
and  1.0- 1.5  mm  groups,  respectively).  This  indicates 
that  NLS  results  are  similar  to  LS  for  small  arteries, 
but  will  significantly  outperform  LS  for  larger  arteries 
(Figs.  7,  8). 

On  the  other  hand,  since  veins  have  very  low  signal 
while  arteries  still  have  high  signal  (provided  compe¬ 
tent  flow  rephasing)  in  SWI  processed  images,  multi¬ 
plying  with  the  SWI  processed  magnitude  mask  may 
also  help  suppress  the  veins  in  the  LS  results  (10,27), 
and  we  have  also  tested  this  method  using  the  same 
data  we  collected  (results  not  shown).  For  large 


Figure  8.  Artery-tissue  CNR  of  different  artery  size  groups. 
Two-tailed  paired  t-test  results  between  LS  and  NLS  are:  P  = 
0.126  for  <1  mm,  P  =  0.395  for  1-1.5  mm,  P  =  0.046  for 
1.5-2  mm  and  P  =  0.035  for  2-2.5  mm. 
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arteries,  masking  with  SWI  images  also  yielded  very 
high  artery-tissue  contrast  and  negative  vein-tissue 
contrast.  Smaller  arteries,  on  the  other  hand,  were 
not  as  equally  enhanced  as  in  LS  or  NLS,  probably 
due  to  the  local  signal  variations  and  increased  noise 
caused  by  phase  masking  process.  Another  potential 
drawback  of  this  SWI  masking  method  would  be  the 
introduction  of  susceptibility  artifacts  commonly  seen 
in  SWI  images  (11),  which  may  suppress  nearby 
arteries  due  to  signal  dropout  especially  at  the  edge  of 
the  brain.  The  SWI  masking  method  may  find  itself 
useful  in  full-dose  CE  MRA  scans,  in  which  both 
arteries  and  veins  display  high  values  in  the  original 
images,  rendering  both  subtraction  methods  less 
effective  in  reducing  venous  signals.  By  integrating 
phase  information,  the  SWI  masking  method  may 
have  the  potential  to  selectively  suppress  the  veins 
even  with  the  presence  of  contrast  agent. 

The  major  drawback  of  the  interleaved  design  of  the 
sequence  is  that  it  still  requires  approximately  the 
same  total  scan  time  as  the  separate  acquisition  of 
both  MRA  and  SWI  data.  However,  compared  to  sepa¬ 
rate  acquisition  of  both  MRA  and  SWI,  our  MRAV 
method  offers  strictly  aligned  images  for  direct  and 
accurate  comparison  between  the  two  vascular  net¬ 
works,  and  a  selectively  enhanced  MRA  map  that  sur¬ 
passes  any  existing  noncontrast-enhanced  MRA  or 
MRAV  methods.  Also  with  proper  adjustment  in  scan¬ 
ning  parameters,  this  method  may  have  the  potential 
for  MRAV  imaging  of  proximal  vessels  such  as  the  ca¬ 
rotid  arteries  and  the  jugular  veins  in  the  neck. 

In  conclusion,  we  have  proposed  an  FR/FD  inter¬ 
leaved  double-echo  MRAV  sequence  which  for  the  first 
time  can  provide  both  fully  flow-compensated  SWI 
results  and  selective  MRA  enhancement  in  a  single 
scan,  while  maintaining  a  good  balance  between  the 
two  contrast  sets.  Three  different  types  of  images,  i.e., 
TOF  MRA,  SWI,  and  flow  dephased,  could  be  obtained 
simultaneously  with  perfect  alignment  to  each  other 
so  that  no  realignment  process  is  necessary.  In  con¬ 
junction  with  NLS,  one  can  obtain  high-resolution 
SWI  and  enhanced  TOF  MRA  without  the  need  for 
contrast  agent,  and  the  process  can  be  easily  auto¬ 
mated  with  minimal  additional  computational  power. 
NLS  offers  a  clean  TOF  angiography  free  of  venous 
contamination  as  if  acquired  with  minimal  saturation 
effects.  For  future  works,  the  acquisition  time  may  be 
further  reduced  to  provide  whole  brain  coverage  using 
a  segmented  echo  planar  imaging  (EPI)  readout 
approach  (28)  or  compress  sensing  (29-31)  since  the 
subtracted  images  contain  highly  sparse  information 
of  arteries.  Also,  a  TONE  pulse  may  be  implemented 
with  proper  flip  angle  range  to  further  improve  the 
enhanced  MRA  results  while  maintaining  high  SWI 
contrast. 
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